/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Daniel Mazin, 8/2004 <mailto:mazin@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MSkyPlot
//
// Create a false source plot. For the Binning in x,y the object MBinning
// "BinningFalseSource" is taken from the paremeter list.
//
// The binning in alpha is currently fixed as 18bins from 0 to 90deg.
//
// If MTime, MObservatory and MPointingPos is available in the paremeter
// list each image is derotated.
//
// MSkyPlot fills a 3D histogram with alpha distribution for
// each (derotated) x and y position.
//
// Only a radius of 1.2deg around the camera center is taken into account.
//
// The displayed histogram is the projection of alpha<fAlphaCut into
// the x,y plain and the projection of alpha>90-fAlphaCut
//
// By using the context menu (find it in a small region between the single
// pads) you can change the AlphaCut 'online'
//
// Each Z-Projection (Alpha-histogram) is scaled such, that the number
// of entries fits the maximum number of entries in all Z-Projections.
// This should correct for the different acceptance in the center
// and at the border of the camera. This, however, produces more noise
// at the border.
//
// Here is a slightly simplified version of the algorithm:
// ------------------------------------------------------
//    MHillas hil; // Taken as second argument in MFillH
//
//    MSrcPosCam src;
//    MHillasSrc hsrc;
//    hsrc.SetSrcPos(&src);
//
//    for (int ix=0; ix<nx; ix++)
//        for (int iy=0; iy<ny; iy++)
//        {
//            TVector2 v(cx[ix], cy[iy]); //cx center of x-bin ix
//            if (rho!=0)                 //cy center of y-bin iy
//                v=v.Rotate(rho);         //image rotation angle
//
//            src.SetXY(v);               //source position in the camera
//
//            if (!hsrc.Calc(hil))        //calc source dependant hillas
//                return;
//
//            //fill absolute alpha into histogram
//            const Double_t alpha = hsrc.GetAlpha();
//            fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
//        }
//    }
//
// Use MHFalse Source like this:
// -----------------------------
//    MFillH fill("MSkyPlot", "MHillas");
// or
//    MSkyPlot hist;
//    hist.SetAlphaCut(12.5);  // The default value
//    hist.SetBgmean(55);      // The default value
//    MFillH fill(&hist, "MHillas");
//
// To be implemented:
// ------------------
//  - a switch to use alpha or |alpha|
//  - taking the binning for alpha from the parlist (binning is
//    currently fixed)
//  - a possibility to change the fit-function (eg using a different TF1)
//  - a possibility to change the fit algorithm (eg which paremeters
//    are fixed at which time)
//  - use a different varaible than alpha
//  - a possiblity to change the algorithm which is used to calculate
//    alpha (or another parameter) is missing (this could be done using
//    a tasklist somewhere...)
//  - a more clever (and faster) algorithm to fill the histogram, eg by
//    calculating alpha once and fill the two coils around the mean
//  - more drawing options...
//  - Move Significance() to a more 'general' place and implement
//    also different algorithms like (Li/Ma)
//  - implement fit for best alpha distribution -- online (Paint)
//  - currently a constant pointing position is assumed in Fill()
//  - the center of rotation need not to be the camera center
//  - ERRORS on each alpha bin to estimate the chi^2 correctly!
//    (sqrt(N)/binwidth) needed for WOlfgangs proposed caluclation
//    of alpha(Li/Ma)
//
//////////////////////////////////////////////////////////////////////////////
#include "MSkyPlot.h"

#include <TF1.h>
#include <TH2.h>
#include <TH1.h>
#include <TGraph.h>
#include <TStyle.h>
#include <TLatex.h>
#include <TCanvas.h>
#include <TPaveText.h>
#include <TStopwatch.h>
#include <TFile.h>

#include "MGeomCam.h"
#include "MSrcPosCam.h"
#include "MReportDrive.h"
#include "MHillasSrc.h"
#include "MHillas.h"
#include "MHillasExt.h"
#include "MNewImagePar.h"
#include "MTime.h"
#include "MObservatory.h"
#include "MPointingPos.h"
#include "MAstroCatalog.h"
#include "MAstroSky2Local.h"
#include "MStarCamTrans.h"
#include "MStatusDisplay.h"
#include "MMath.h"
#include "MSupercuts.h"

#include "MBinning.h"
#include "MParList.h"

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

ClassImp(MSkyPlot);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor
//
MSkyPlot::MSkyPlot(const char *name, const char *title)
    : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1) 
//    fAlphaCut(12.5), BgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1)
{

*fLog << warn << "entering default constructor in MSkyPlot" << endl; 
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MSkyPlot";
    fTitle = title ? title : "sky plot vs x, y";

    fSetCenter=kTRUE;

    fHistSignif.SetDirectory(NULL);
    fHistNexcess.SetDirectory(NULL);
    fHistOn.SetDirectory(NULL);

    fHistSignif.SetName("SkyPlotSignif");
    fHistSignif.SetTitle("Sky Plot of significance vs x, y");
    fHistSignif.SetXTitle("x [\\circ]");
    fHistSignif.SetYTitle("y [\\circ]");

    fHistNexcess.SetName("SkyPlotNexcess");
    fHistNexcess.SetTitle("Sky Plot of number of excess vs x, y");
    fHistNexcess.SetXTitle("x [\\circ]");
    fHistNexcess.SetYTitle("y [\\circ]");

    fHistOn.SetName("SkyPlotOn");
    fHistOn.SetTitle("Sky Plot of number of On events vs x, y");
    fHistOn.SetXTitle("x [\\circ]");
    fHistOn.SetYTitle("y [\\circ]");

    // set some values for the sky plot geometry:
    fMinXGrid    =-1.;  		// [deg] 
    fMaxXGrid    = 1.;             //  [deg] , right edge of the skyplot
    fMinYGrid    =-1.;            //  [deg] , upper edge of the skyplot
    fMaxYGrid    = 1.;             //  [deg] , lower edge of the skyplot
    fBinStepGrid = 0.05;         //  [deg], 
    fAlphaONMax  = 5.;           //  [deg] , upper cut for alpha ON region in the alpha plot
                      // [deg], ON region in the alpha plot, maybe 5 deg is better
                      // NOTE: up to now only values of 5, 10, 15, 20 degrees are possible
// ra,dec lines from wolfgang:
    fGridBinning = 0.50;   // degrees
    fGridFineBin = 0.01;   // degrees


    // some filter cuts:
    fSizeMin    = 2000.;             // min size in photons
    fSizeMax    = 100000.;             // max size in photons
    fLeakMax    = 0.25;             // leakmax  in per cent
    fMaxDist = 1.4;          // dist max cut (ever possible)
    fMinDist = 0.1;          // dist max cut (ever possible)

    fNumBinsAlpha   = 36; 	 // number of bins for alpha histograms
    fAlphaLeftEdge  = 0.;        // [deg] left edge
    fAlphaRightEdge = 90.;       // [deg] left edge

    fAlphaBgLow     = 30.;
    fAlphaBgUp      = 90.;

    {
// SET DEFAULT VALUES HERE	
    fLengthUp[0] =  0.2;
    fLengthUp[1] =  0.0;
    fLengthUp[2] =  0.0;
    fLengthUp[3] =  0.0;
    fLengthUp[4] =  0.0;
    fLengthUp[5] =  0.0;
    fLengthUp[6] =  0.0;
    fLengthUp[7] =  0.0;

    fLengthLo[0] =  0.;
    fLengthLo[1] =  0.;
    fLengthLo[2] =  0.;
    fLengthLo[3] =  0.;
    fLengthLo[4] =  0.;
    fLengthLo[5] =  0.;
    fLengthLo[6] =  0.;
    fLengthLo[7] =  0.;

    fWidthUp[0] =  0.1;
    fWidthUp[1] =  0.0;
    fWidthUp[2] =  0.0;
    fWidthUp[3] =  0.0;
    fWidthUp[4] =  0.0;
    fWidthUp[5] =  0.0;
    fWidthUp[6] =  0.0;
    fWidthUp[7] =  0.0;

    fWidthLo[0] =  0.;
    fWidthLo[1] =  0.;
    fWidthLo[2] =  0.;
    fWidthLo[3] =  0.;
    fWidthLo[4] =  0.;
    fWidthLo[5] =  0.;
    fWidthLo[6] =  0.;
    fWidthLo[7] =  0.;

    fDistUp[0] =  1.e10;
    fDistUp[1] =  0.0;
    fDistUp[2] =  0.0;
    fDistUp[3] =  0.0;
    fDistUp[4] =  0.0;
    fDistUp[5] =  0.0;
    fDistUp[6] =  0.0;
    fDistUp[7] =  0.0;

    fDistLo[0] =  0.0;
    fDistLo[1] =  0.0;
    fDistLo[2] =  0.0;
    fDistLo[3] =  0.0;
    fDistLo[4] =  0.0;
    fDistLo[5] =  0.0;
    fDistLo[6] =  0.0;
    fDistLo[7] =  0.0;
    }
}

void MSkyPlot::ReadCuts(const TString parSCinit="OptSCParametersONOFFThetaRange0_1570mRad.root")
{
    //--------------------------------
    // create container for the supercut parameters
    // and set them to their initial values
    MSupercuts super;

  // read the cuts coefficients
    TFile inparam(parSCinit);
    super.Read("MSupercuts");
    inparam.Close();
    *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from file "
          << parSCinit << endl;

    TArrayD params = super.GetParameters();
    TArrayD steps  = super.GetStepsizes();
 // TMP2
    if (params.GetSize() == steps.GetSize())
    {
            *fLog << "SC parameters and Setps are: " << endl;
            for (Int_t z = 0; z < params.GetSize(); z++)
            {
                cout << params[z] << setw(20) << steps[z] << endl;
            }
    }
 // ENDTMP2
    for (Int_t i=0; i<8; i++)
    {
	fLengthUp[i]  = params[i];
        fLengthLo[i] = params[i+8];
	fWidthUp[i]   = params[i+16];
	fWidthLo[i]  = params[i+24];
        fDistUp[i]    = params[i+32];
        fDistLo[i]   = params[i+40];
    }
}

void MSkyPlot::SetSkyPlot(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t step)
{
    Float_t temp;

    if (xmax<xmin)
    {
        *fLog << warn << "SetSkyPlot: xmax is smaller xmin ... exchanging them." << endl;
	temp = xmax;
	xmax = xmin;
	xmin = temp;
    }	

    if (ymax<ymin)
    {
        *fLog << warn << "SetSkyPlot: ymax is smaller ymin ... exchanging them." << endl;
	temp = ymax;
	ymax = ymin;
	ymin = temp;
    }	
    
    if (step<0)
        *fLog << warn << "SetSkyPlot: step<0... taking absolute value." << endl;

    fBinStepGrid = TMath::Abs(step);
    fMinXGrid    = xmin;
    fMaxXGrid    = xmax;
    fMinYGrid    = ymin;
    fMaxYGrid    = ymax;
}

// --------------------------------------------------------------------------
//
// Set the alpha cut (|alpha|<fAlphaCut) which is use to estimate the
// number of excess events
//
void MSkyPlot::SetAlphaCut(Float_t alpha)
{
    if (alpha<0)
        *fLog << warn << "Alpha<0... taking absolute value." << endl;

    fAlphaONMax = TMath::Abs(alpha);
}

// calculate the values for the cuts:
Double_t MSkyPlot::CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2)
{

    Double_t  limit = a[0] + a[1] * dd2 +
                      ls  * (a[3] + a[4] * dd2) +
                      ls2 * (a[6] + a[7] * dd2);

     return limit;
}



// --------------------------------------------------------------------------
//
// Set binnings (takes BinningFalseSource) and prepare filling of the
// histogram.
//
// Also search for MTime, MObservatory and MPointingPos
// 
Int_t MSkyPlot::PreProcess(MParList *plist)
{
*fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl; 
    fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
    if (!fGeomCam)
    {
        *fLog << err << "MGeomCam not found... aborting." << endl;
        return kFALSE;
    }

    fMm2Deg = fGeomCam->GetConvMm2Deg();

    fRepDrive = (MReportDrive*)plist->FindObject(AddSerialNumber("MReportDrive"));
    if (!fRepDrive)
        *fLog << warn << "MReportDrive not found... no sky map." << endl;


    fSrcPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
    if (!fSrcPosCam)
        *fLog << warn << "MSrcPosCam not found... no sky map." << endl;

    fPntPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MPntPosCam"));
    if (!fPntPosCam)
        *fLog << warn << "MPntPosCam not found... no sky map." << endl;

    fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
    if (!fPointPos)
        *fLog << warn << "MPointingPos not found... no sky map." << endl;

    fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime"));
    if (!fTime)
        *fLog << warn << "MTime not found... no sky map." << endl;

    fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
    if (!fObservatory)
        *fLog << warn << "MObservatory not found... no sky map." << endl;


    fHillas = (MHillas*)plist->FindObject(AddSerialNumber("MHillas"));
    if (!fHillas)
        *fLog << err << "MHillas not found... no sky map." << endl;

    fHillasExt = (MHillasExt*)plist->FindObject(AddSerialNumber("MHillasExt"));
    if (!fHillasExt)
        *fLog << err << "MHillasExt not found... no sky map." << endl;

    fHillasSrc = (MHillasSrc*)plist->FindObject(AddSerialNumber("MHillasSrc"));
    if (!fHillasSrc)
        *fLog << err << "MHillasSrc not found... no sky map." << endl;

    fNewImagePar = (MNewImagePar*)plist->FindObject(AddSerialNumber("MNewImagePar"));
    if (!fNewImagePar)
        *fLog << err << "MNewImagePar not found... no sky map." << endl;



    // FIXME: Because the pointing position could change we must check
    // for the current pointing position and add a offset in the
    // Fill function!
    kSaveAlphaPlots=kTRUE;
    kSaveSkyPlots=kTRUE;

   // prepare skyplot
    fNumStepsX     =  (int) ((fMaxXGrid - fMinXGrid) / fBinStepGrid + 1.5);
    fNumStepsY     =  (int) ((fMaxYGrid - fMinYGrid) / fBinStepGrid + 1.5);
    fNumalphahist  =  (int) (fNumStepsX * fNumStepsY  + 0.5);

    fHistSignif.SetName("SPSignif2ndOrder");
    fHistSignif.SetTitle("Sky Plot of significance (2nd order fit)");
    fHistSignif.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
    fHistSignif.SetDirectory(NULL);
    fHistSignif.SetFillStyle(4000);
    fHistSignif.UseCurrentStyle();

    fHistNexcess.SetName("SPNex2ndOrder");
    fHistNexcess.SetTitle("Sky Plot of Number of excess events (2nd order fit)");
    fHistNexcess.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
    fHistNexcess.SetDirectory(NULL);
    fHistNexcess.SetFillStyle(4000);
    fHistNexcess.UseCurrentStyle();

    fHistOn.SetName("SPOnCounted");
    fHistOn.SetTitle("Sky Plot of ON events (counted)");
    fHistOn.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
    fHistOn.SetDirectory(NULL);
    fHistOn.SetFillStyle(4000);
    fHistOn.UseCurrentStyle();


    // prepare alpha plots
    	// temp histogram for an alpha plot
    TH1I *temp = new TH1I("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
    temp->SetDirectory(NULL);
    fHistAlphaBinWidth = temp->GetBinWidth(1);
    // create alpha histograms
    for (Int_t i=0; i< fNumalphahist; i++) {
    	fHistAlpha.push_back(*temp);
	fHistAlpha[i].SetDirectory(NULL);
    }

    delete temp;
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the histogram. For details see the code or the class description
// 
Int_t MSkyPlot::Process()
{
// check whether MPointingPos comtains something:
  if (TMath::Abs(fPointPos->GetRa())<1e-3 && TMath::Abs(fPointPos->GetDec())<1e-3) 
  {
     *fLog << warn << "MPointingPos is not filled ... event skipped" << endl;
     return kTRUE;
  }

    // Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside)
  if (fSetCenter==kTRUE)
  {
     fRa0  = fRepDrive->GetRa();  // [h]
     fDec0 = fRepDrive->GetDec(); // [deg]
     fSetCenter=kFALSE;
  }


// some filter cuts:
  if ( fHillas->GetSize() > fSizeMin && fHillas->GetSize() < fSizeMax && fNewImagePar->GetLeakage1() < fLeakMax)
  {

    Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam;
    Double_t dx, dy, beta, tanbeta, alphapar, distpar;
    Double_t logsize, lgsize, lgsize2, dist2; 
    const Double_t log3000 = log(3000.);
    const Float_t fDistOffset = 0.9;

    // Get camera rotation angle
    Double_t rho = 0;
    if (fTime && fObservatory && fPointPos)
{
        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
*fLog << inf << " rho = " << rho << ", Zd = " << fPointPos->GetZd() << 
                ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
}
    //if (fPointPos)
    //    rho = fPointPos->RotationAngle(*fObservatory);
/*
//TEMP
//        theta = mcevt->GetTelescopeTheta();
    Double_t theta, phi, sing, cosg;
    theta = fPointPos->GetZd()*TMath::Pi()/180.;
    phi = fPointPos->GetAz()*TMath::Pi()/180.;
//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
    fObservatory->RotationAngle(theta, phi, sing, cosg);
      
// conclusion: diffference in rho = 7 deg
  *fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ",  "
        << phi << ",  " << cosg << ",  " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl;

    Double_t a1 =  0.876627;
    Double_t a3 = -0.481171;
    Double_t thetaTel=theta, phiTel=phi;

    Double_t denom =  1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
                              ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) *
                              ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) )   );
    Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
    Double_t sinal =    a1 * sin(phiTel) * denom;

  *fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ",  "
        << phiTel << ",  " << cosal << ",  " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl;



// END TEMP
*/

    // => coord. system: xtilde, ytilde with center in RA_0, DEC_0 

// Get Rot. Angle:
	cosgam = TMath::Cos(rho);
	singam = TMath::Sin(rho);
// Get x_0, y_0 for RA_0,DEC_0 of the current event
        x_0 = fPntPosCam->GetX()*fMm2Deg; 
        y_0 = fPntPosCam->GetY()*fMm2Deg; 

    Int_t index = 0;  // index for alpha histograms
    // loop over xtilde
    for (Int_t gridy = 0; gridy < fNumStepsY; gridy++)   
    {
	ysource = fMinYGrid + gridy * fBinStepGrid;
        // loop over ytilde
	for (Int_t gridx = 0; gridx < fNumStepsX; gridx++)
	{
	   xsource = fMinXGrid + gridx * fBinStepGrid;
 /*     derotation    : rotate  into camera coordinates */
 /*     formel: (x_cam)    (cos(gam)  -sin(gam))   (xtilde)   (x_0)
                (     )  = (                   ) * (      ) + (   ) 
	        (y_cam)    (sin(gam)   sin(gam))   (ytilde)   (y_0)
*/
           xsourcam = cosgam * xsource - singam * ysource + x_0;
           ysourcam = singam * xsource + cosgam * ysource + y_0;
 /*    end derotatiom    */

       /* calculate ALPHA und DIST according to the source position */
           dx = fHillas->GetMeanX()*fMm2Deg - xsourcam;
           dy = fHillas->GetMeanY()*fMm2Deg - ysourcam;
           tanbeta = dy / dx ;
           beta = TMath::ATan(tanbeta);
           alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi();
           distpar = sqrt( dy*dy + dx*dx );
           if(alphapar >  90.) alphapar -= 180.;
           if(alphapar < -90.) alphapar += 180.;
           alphapar = TMath::Abs(alphapar);

// apply cuts
           logsize = log(fHillas->GetSize());
           lgsize = logsize-log3000;
           lgsize2 = lgsize*lgsize;
           dist2 = (distpar-fDistOffset)*(distpar-fDistOffset);

           if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
                (fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2))
           if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) &&
                (fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2))
           if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) &&
                      distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) &&
                      distpar < fMaxDist && distpar > fMinDist)
	   {
// gamma candidates!
		fHistAlpha[index].Fill(alphapar);
		index++;
	   }
	}
     }
  }

    return kTRUE;
}

// calculate number of events below alphacut, number of excess events, significance

Int_t MSkyPlot::PostProcess()
{

    Int_t nrow, ncolumn;
    Double_t Non, chisquarefit, numdegfreed, Noff_fit, Nex, Sign; 
    const Int_t onbinsalpha = TMath::Nint(fAlphaONMax/fHistAlphaBinWidth);

// fit function for the background
    TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", fAlphaBgLow, fAlphaBgUp);
    fitbgpar->SetLineColor(2);

// number degrees of freedom:
    numdegfreed = (fAlphaBgUp - fAlphaBgLow) / fHistAlphaBinWidth - 2.; 

    vector <TH1I>::iterator alpha_iterator;
    int index2 = 0;  // index of the TH2F histograms

    alpha_iterator = fHistAlpha.begin();

    while( alpha_iterator != fHistAlpha.end() ) {
	 // find the bin in the 2dim histogram
         nrow    = index2/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
         ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)

         // number of ON events below fAlphaONMax
         Non = 0.;
         for(Int_t i=1; i<=onbinsalpha;i++) Non += (*alpha_iterator).GetBinContent(i);

         fHistOn.SetBinContent(ncolumn, nrow, Non);		// fill in the fHistOn

         // fit parabola to a background region
         (*alpha_iterator).Fit(fitbgpar,"NWR");  // NWR OK?????????????????????????
         // get Chi2
         chisquarefit = fitbgpar->GetChisquare();
         if (chisquarefit/numdegfreed<2. && chisquarefit/numdegfreed>0.01);
         else  *fLog << warn << "Fit is bad, chi2/ndf = " << chisquarefit/numdegfreed << endl;

         // estimate Noff from the fit:
       	 Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * pow(fAlphaONMax,3.) +
       	            (fitbgpar->GetParameter(1)) * fAlphaONMax) /  fHistAlphaBinWidth;

	 Nex = Non - Noff_fit;

         fHistNexcess.SetBinContent(ncolumn, nrow, Nex);	// fill in the fHistNexcess

         if (Noff_fit<0.) Sign = 0.;
//         else Sign = LiMa17(Non,Noff_fit,1.);
         else Sign = MMath::SignificanceLiMaSigned(Non, Noff_fit, 1.);

         fHistSignif.SetBinContent(ncolumn, nrow, Sign);	// fill in the fHistSignif

         alpha_iterator++;
         index2++;
    }

//temp
   Double_t decl, hang;
   MStarCamTrans mstarc(*fGeomCam, *fObservatory);
   mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang);

*fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl;
*fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl;
//endtemp

// save alpha plots:
  if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots();

// save sky plots:
  if(kSaveSkyPlots==kTRUE) SaveSkyPlots();

  fHistAlpha.clear();
  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude
// is set to 9, while the fov is equal to the current fov of the false
// source plot.
//
TObject *MSkyPlot::GetCatalog()
{
    const Double_t maxr = 0.98*TMath::Abs(fHistSignif.GetBinCenter(1));

    // Create catalog...
    MAstroCatalog stars;
    stars.SetLimMag(9);
    stars.SetGuiActive(kFALSE);
    stars.SetRadiusFOV(maxr);
    stars.SetRaDec(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad());
    stars.ReadBSC("bsc5.dat");

    TObject *o = (MAstroCatalog*)stars.Clone();
    o->SetBit(kCanDelete);

    return o;
}


void MSkyPlot::SaveSkyPlots(const TString skyplotfilename="skyplots.root")
{
   TFile rootfile(skyplotfilename, "RECREATE",
                   "sky plots in some variations");
     fHistSignif.Write();
     fHistNexcess.Write();
     fHistOn.Write();

// from Wolfgang: 
   //--------------------------------------------------------------------
    // Dec-Hour lines
    Double_t rho, sinrho, cosrho;
//    if (fTime && fObservatory && fPointPos)
        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
    sinrho=TMath::Sin(rho);
    cosrho=TMath::Cos(rho);

    Double_t theta, phi, sing, cosg;
    theta = fPointPos->GetZd()*TMath::Pi()/180.;
    phi = fPointPos->GetAz()*TMath::Pi()/180.;
//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
    fObservatory->RotationAngle(theta, phi, sing, cosg);

*fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl;
*fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl;

    Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0;     //  [mm]
    Double_t gridbinning = fGridBinning;
    Double_t gridfinebin = fGridFineBin;
    Int_t    numgridlines = (Int_t)(4.0/gridbinning);
    Double_t aberr = 1.07;
    Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;

    Double_t declin, hangle;  // [deg], [h]
    MStarCamTrans mstarc(*fGeomCam, *fObservatory);
    mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);

    TGraph *graph1;
    TLatex *pix;

    Char_t tit[100];
    Double_t xtxt;
    Double_t ytxt;

    Double_t xlo = -534.0 * mmtodeg;
    Double_t xup =  534.0 * mmtodeg;

    Double_t ylo = -534.0 * mmtodeg;
    Double_t yup =  534.0 * mmtodeg;

    Double_t xx, yy;


    // direction for the center of the camera
    Double_t dec0 = declin;
    Double_t h0   = hangle*360./24.;
    //trans.LocToCel(theta0, phi0, dec0, h0);

    TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600);
    c1->cd(1);
    fHistOn.Draw("colz");

    //-----   lines for fixed dec   ------------------------------------

    const Int_t Ndec = numgridlines;
    Double_t dec[Ndec];
    Double_t ddec = gridbinning;

    Int_t nh = (Int_t)(4.0/gridfinebin);
    const Int_t Nh   = nh+1;
    Double_t h[Nh];
    Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
    if ( dh > 360.0/((Double_t)(Nh-1)) )
      dh = 360.0/((Double_t)(Nh-1));


    for (Int_t j=0; j<Ndec; j++)
    {
      dec[j] = dec0 + ((Double_t)j)*ddec
                        - ((Double_t)(Ndec/2)-1.0)*ddec;
    }

    for (Int_t k=0; k<Nh; k++)
    {
      h[k] = h0 + ((Double_t)k)*dh - ((Double_t)(Nh/2)-1.0)*dh;
    }

    Double_t xh[Nh];
    Double_t yh[Nh];



    for (Int_t j=0; j<Ndec; j++)
    {
      if (fabs(dec[j]) > 90.0) continue;

      for (Int_t k=0; k<Nh; k++)
      {
        Double_t hh0 = h0   *24.0/360.0;
        Double_t hx = h[k]*24.0/360.0;
        mstarc.Cel0CelToCam(dec0, hh0, dec[j], hx,
                           xx, yy);
        xx = xx * mmtodeg * aberr;
        yy = yy * mmtodeg * aberr;
//        xh[k] = xx * mmtodeg * aberr;
//        yh[k] = yy * mmtodeg * aberr;
        xh[k] = cosg * xx + sing * yy;
        yh[k] =-sing * xx + cosg * yy;
//        xh[k] = cosrho * xx + sinrho * yy;
//        yh[k] =-sinrho * xx + cosrho * yy;


        //gLog << "dec0, h0 = " << dec0 << ",  " << h0
        //     << " : dec, h, xh, yh = " << dec[j] << ",  "
        //     << h[k] << ";   " << xh[k] << ",  " << yh[k] << endl;
      }

//      c1->cd(2);
      graph1 = new TGraph(Nh,xh,yh);
      graph1->SetLineColor(j+1);
      graph1->SetLineStyle(1);
      graph1->SetMarkerColor(j+1);
      graph1->SetMarkerSize(.2);
      graph1->SetMarkerStyle(20);
      graph1->Draw("PL");
      //delete graph1;

      sprintf(tit,"Dec = %6.2f", dec[j]);
      xtxt = xlo + (xup-xlo)*0.1;
      ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0;
      pix = new TLatex(xtxt, ytxt, tit);
      pix->SetTextColor(j+1);
      pix->SetTextSize(.03);
      pix->Draw("");
      //delete pix;

    }

    //-----   lines for fixed hour angle   ------------------------------------

    Int_t ndec1 = (Int_t)(4.0/gridfinebin);
    const Int_t Ndec1 = ndec1;
    Double_t dec1[Ndec1];
    Double_t ddec1 = gridfinebin;

    const Int_t Nh1   = numgridlines;
    Double_t h1[Nh1];
    Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
    if ( dh1 > 360.0/((Double_t)Nh1) )
      dh1 = 360.0/((Double_t)Nh1);


    for (Int_t j=0; j<Ndec1; j++)
    {
      dec1[j] = dec0 + ((Double_t)j)*ddec1
                        - ((Double_t)(Ndec1/2)-1.0)*ddec1;
    }

    for (Int_t k=0; k<Nh1; k++)
    {
      h1[k] = h0 + ((Double_t)k)*dh1 - ((Double_t)(Nh1/2)-1.0)*dh1;
    }


    Double_t xd[Ndec1];
    Double_t yd[Ndec1];

    for (Int_t k=0; k<Nh1; k++)
    {
      Int_t count = 0;
      for (Int_t j=0; j<Ndec1; j++)
      {
        if (fabs(dec1[j]) > 90.0) continue;

        Double_t hh0 = h0   *24.0/360.0;
        Double_t hhx = h1[k]*24.0/360.0;
        mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx,
                           xx, yy);
//        xd[count] = xx * mmtodeg * aberr;
//        yd[count] = yy * mmtodeg * aberr;

        xx = xx * mmtodeg * aberr;
        yy = yy * mmtodeg * aberr;
        xd[count] = cosg * xx + sing * yy;
        yd[count] =-sing * xx + cosg * yy;

        //gLog << "dec0, h0 = " << dec0 << ",  " << h0
        //     << " : dec1, h1, xd, yd = " << dec1[j] << ",  "
        //     << h1[k] << ";   " << xd[count] << ",  " << yd[count] << endl;

        count++;
      }

//      c1->cd(2);
      graph1 = new TGraph(count,xd,yd);
      graph1->SetLineColor(k+1);
      graph1->SetLineStyle(2);
      graph1->SetMarkerColor(k+1);
      graph1->SetMarkerSize(.2);
      graph1->SetMarkerStyle(20);
      graph1->Draw("PL");
      //delete graph1;

      sprintf(tit,"H = %6.2f", h1[k]);
      Double_t xtxt = xlo + (xup-xlo)*0.75;
      Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0;
      pix = new TLatex(xtxt, ytxt, tit);
      pix->SetTextColor(k+1);
      pix->SetTextSize(.03);
      pix->Draw("");
      //delete pix;

    }

//    c1->cd(2);
    sprintf(tit,"Dec0 = %6.2f    H0 = %6.2f", dec0, h0);
    xtxt = xlo + (xup-xlo)*0.05;
    ytxt = ylo + (yup-ylo)*0.92;
    pix = new TLatex(xtxt, ytxt, tit);
    pix->SetTextColor(1);
    pix->SetTextSize(.06);
    pix->Draw("");
    //delete pix;

    sprintf(tit,"              [deg]             [deg]");
    xtxt = xlo + (xup-xlo)*0.05;
    ytxt = ylo + (yup-ylo)*0.86;
    pix = new TLatex(xtxt, ytxt, tit);
    pix->SetTextColor(1);
    pix->SetTextSize(.06);
    pix->Draw("");

    c1->Write();

    TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600);
    c2->Divide(1,2);
    c2->SetBorderMode(0);
    c2->cd(1);
    fHistSignif.Draw("colz");
    c2->cd(2);
    fHistNexcess.Draw("colz");
    c2->Write();

   rootfile.Close();
   delete c1;
   delete c2;
   delete graph1;
   delete pix;
}

void MSkyPlot::SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root")
{
   TFile rootfile(alphaplotfilename, "RECREATE",
                   "all the alpha plots");

    vector <TH1I>::iterator alpha_iterator;
    int index = 0;  // index of the TH2F histograms
    Char_t strtitle[100];
    Char_t strname[100];
    Float_t xpos, ypos, signif, nex;
    Int_t nrow, ncolumn, non;


    alpha_iterator = fHistAlpha.begin();

    while( alpha_iterator != fHistAlpha.end() ) {

          nrow    = index/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
          ncolumn = TMath::Nint(fmod(index,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
          xpos    = fMinXGrid + fBinStepGrid*(ncolumn-1);
          ypos    = fMinYGrid + fBinStepGrid*(nrow-1);
          non     = TMath::Nint(fHistOn.GetBinContent(ncolumn,nrow));
          nex     = fHistNexcess.GetBinContent(ncolumn,nrow);
          signif  = fHistSignif.GetBinContent(ncolumn,nrow);

          sprintf(strname,"AlphaPlotX%.2fY%.2f", xpos, ypos); 
          sprintf(strtitle,"for x= %.2f deg, y= %.2f deg: S= %.2f sigma, Nex= %.2f, Non= %d", 
			xpos, ypos, signif, nex, non);
           
	  (*alpha_iterator).SetName(strname);
	  (*alpha_iterator).SetTitle(strtitle);
          (*alpha_iterator).Write();
	  alpha_iterator++;
          index++;
    }

   rootfile.Close();
}

// --------------------------------------------------------------------------
//
// Draw the histogram
//
void MSkyPlot::Draw(Option_t *opt)
{
/*
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);

    AppendPad("");

    pad->Divide(1, 2, 0, 0.03);

    TObject *catalog = GetCatalog();

    // Initialize upper part
    pad->cd(1);
    gPad->SetBorderMode(0);
    gPad->Divide(3, 1);
*/
/*
    // PAD #1
    pad->GetPad(1)->cd(1);
    gPad->SetBorderMode(0);
    fHist.GetZaxis()->SetRangeUser(0,fAlphaONMax);
    TH1 *h3 = fHist.Project3D("yx_on");
    fHist.GetZaxis()->SetRange(0,9999);
    h3->SetDirectory(NULL);
    h3->SetXTitle(fHist.GetXaxis()->GetTitle());
    h3->SetYTitle(fHist.GetYaxis()->GetTitle());
    h3->Draw("colz");
    h3->SetBit(kCanDelete);
    catalog->Draw("mirror same");

    // PAD #2
    pad->GetPad(1)->cd(2);
    gPad->SetBorderMode(0);
    fHist.GetZaxis()->SetRange(0,0);
    TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning....
    fHist.GetZaxis()->SetRange(0,9999);
    h4->SetTitle("Significance");
    h4->SetDirectory(NULL);
    h4->SetXTitle(fHist.GetXaxis()->GetTitle());
    h4->SetYTitle(fHist.GetYaxis()->GetTitle());
    h4->Draw("colz");
    h4->SetBit(kCanDelete);
    catalog->Draw("mirror same");

    // PAD #3
    pad->GetPad(1)->cd(3);
    gPad->SetBorderMode(0);
    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
    TH1 *h2 = fHist.Project3D("yx_off");
    h2->SetDirectory(NULL);
    h2->SetXTitle(fHist.GetXaxis()->GetTitle());
    h2->SetYTitle(fHist.GetYaxis()->GetTitle());
    h2->Draw("colz");
    h2->SetBit(kCanDelete);
    catalog->Draw("mirror same");


    // Initialize lower part
    pad->cd(2);
    gPad->SetBorderMode(0);
    gPad->Divide(3, 1);

    // PAD #4
    pad->GetPad(2)->cd(1);
    gPad->SetBorderMode(0);
    TH1 *h1 = fHist.ProjectionZ("Alpha");
    h1->SetDirectory(NULL);
    h1->SetTitle("Distribution of \\alpha");
    h1->SetXTitle(fHist.GetZaxis()->GetTitle());
    h1->SetYTitle("Counts");
    h1->Draw(opt);
    h1->SetBit(kCanDelete);

    // PAD #5
    pad->GetPad(2)->cd(2);
    gPad->SetBorderMode(0);
    TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff");
    h5->Add(h2, -1);
    h5->SetTitle("Difference of on- and off-distribution");
    h5->SetDirectory(NULL);
    h5->Draw("colz");
    h5->SetBit(kCanDelete);
    catalog->Draw("mirror same");

    // PAD #6
    pad->GetPad(2)->cd(3);
    gPad->SetBorderMode(0);
    TH1 *h0 = fHist.Project3D("yx_all");
    h0->SetDirectory(NULL);
    h0->SetXTitle(fHist.GetXaxis()->GetTitle());
    h0->SetYTitle(fHist.GetYaxis()->GetTitle());
    h0->Draw("colz");
    h0->SetBit(kCanDelete);
    catalog->Draw("mirror same");
*/
}

// --------------------------------------------------------------------------
//
// Everything which is in the main pad belongs to this class!
//
Int_t MSkyPlot::DistancetoPrimitive(Int_t px, Int_t py)
{
    return 0;
}

// --------------------------------------------------------------------------
//
// Set all sub-pads to Modified()
//
/*
void MSkyPlot::Modified()
{
    if (!gPad)
        return;

    TVirtualPad *padsave = gPad;
    padsave->Modified();
    padsave->GetPad(1)->cd(1);
    gPad->Modified();
    padsave->GetPad(1)->cd(3);
    gPad->Modified();
    padsave->GetPad(2)->cd(1);
    gPad->Modified();
    padsave->GetPad(2)->cd(2);
    gPad->Modified();
    padsave->GetPad(2)->cd(3);
    gPad->Modified();
    gPad->cd();
}
*/
// --------------------------------------------------------------------------
//
// This is a preliminary implementation of a alpha-fit procedure for
// all possible source positions. It will be moved into its own
// more powerfull class soon.
//
// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
//   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
// or
//   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
//
// Parameter [1] is fixed to 0 while the alpha peak should be
// symmetric around alpha=0.
//
// Parameter [4] is fixed to 0 because the first derivative at
// alpha=0 should be 0, too.
//
// In a first step the background is fitted between bgmin and bgmax,
// while the parameters [0]=0 and [2]=1 are fixed.
//
// In a second step the signal region (alpha<sigmax) is fittet using
// the whole function with parameters [1], [3], [4] and [5] fixed.
//
// The number of excess and background events are calculated as
//   s = int(0, sigint, gaus(0)+pol2(3))
//   b = int(0, sigint,         pol2(3))
//
// The Significance is calculated using the Significance() member
// function.
//

// I might need this
/*
    TCanvas *c=new TCanvas;

    gStyle->SetPalette(1, 0);

    c->Divide(3,2, 0, 0);
    c->cd(1);
    gPad->SetBorderMode(0);
    hists->Draw("colz");
    hists->SetBit(kCanDelete);
    catalog->Draw("mirror same");
    c->cd(2);
    gPad->SetBorderMode(0);
    hist->Draw("colz");
    hist->SetBit(kCanDelete);
    catalog->Draw("mirror same");
    c->cd(3);
    gPad->SetBorderMode(0);
    histb->Draw("colz");
    histb->SetBit(kCanDelete);
    catalog->Draw("mirror same");
    c->cd(4);
    gPad->Divide(1,3, 0, 0);
    TVirtualPad *p=gPad;
    p->SetBorderMode(0);
    p->cd(1);
    gPad->SetBorderMode(0);
    h0b.DrawCopy();
    h0a.DrawCopy("same");
    p->cd(2);
    gPad->SetBorderMode(0);
    h3.DrawCopy();
    p->cd(3);
    gPad->SetBorderMode(0);
    h2.DrawCopy();

*/