/* ======================================================================== *\
!
! *
! * 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): Marcos Lopez 03/2004 <mailto:marcos@gae.ucm.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//   MExtrapolatePointingPos
//
//   In the PreProcess, read the drive report and store it in an TSpline.
//   In the Process, use the TSpline to calculate the PointingPos for the 
//   time of the current event.
// 
//  Input Containers:
//   MRawEvtData
//
//  Output Containers:
//   MPointingPos
//
//
/////////////////////////////////////////////////////////////////////////////

#include "MExtrapolatePointingPos.h"

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

#include <TString.h>
#include <TSpline.h>

#include "MTaskList.h"
#include "MParList.h"
#include "MEvtLoop.h"
#include "MReadReports.h"
#include "MReportDrive.h"
#include "MPointingPos.h"
#include "MTime.h"
#include "MRawRunHeader.h"

#include <TCanvas.h>

ClassImp(MExtrapolatePointingPos);

using namespace std;


// ---------------------------------------------------------------------------
//
//
//
void MExtrapolatePointingPos::ReadDriveReport(const TString filename)
{

    //
    // ParList
    // -------
    //
    MParList  plist;
    MTaskList tlist;
    plist.AddToList(&tlist);


    //
    // TaskList
    // --------
    //
    // Reads the trees of the root file and the analysed branches
    MReadReports read;
    read.AddTree("Drive");
    read.AddFile(filename);     // after the reading of the trees!!!
    read.AddToBranchList("MReportDrive.*");
    
    // Add all the task to the task list
    tlist.AddToList(&read);
 
    
    
    //
    // EventLoop
    // ---------
    //
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);

    
    //
    // Execute your analysis
    //
    if (!evtloop.PreProcess())
        return;


    //
    // Store the ReportTime, CurrentZd, CurrentAz, Ra, Dec and Ha into arrays
    //
    TArrayD ReportTime(10000);
    TArrayD CurrentZd(10000);
    TArrayD CurrentAz(10000); 
    TArrayD NominalZd(10000);
    TArrayD NominalAz(10000);
    TArrayD Ra(10000);
    TArrayD Dec(10000);
    TArrayD Ha(10000);
    
    
    Int_t n=0;

    while (tlist.Process())
    { 
	
	MReportDrive* report = (MReportDrive*)plist.FindObject("MReportDrive");
	MTime* reporttime = (MTime*)plist.FindObject("MTimeDrive");

	//
	// Sometimes there are two reports with the same time
	//
	if (reporttime->GetTime() == ReportTime[n-1])
	{ 
	    //cout << " *********** Error " << endl;
	    continue;
	}

	ReportTime[n] = reporttime->GetTime();
	CurrentZd[n] = report->GetCurrentZd();
	CurrentAz[n] = report->GetCurrentAz();
	NominalZd[n] = report->GetNominalZd();
	NominalAz[n] = report->GetNominalAz();
	Ra[n] = report->GetRa();
	Dec[n] = report->GetDec();
	Ha[n] = report->GetHa();

	n++;
    }

    //tlist.PrintStatistics();


    //
    // Update the number of entries
    //
    ReportTime.Set(n);
    CurrentZd.Set(n);
    CurrentAz.Set(n); 
    NominalZd.Set(n);
    NominalAz.Set(n);
    Ra.Set(n);
    Dec.Set(n);
    Ha.Set(n);
    


    //   for(int i=0;i<time.GetSize();i++)
//       {
//  	 cout << i << " " << time[i] << " " << Zd[i] << endl;
//       }
     

     

    fSplineZd = new TSpline3("zenith",
			     ReportTime.GetArray(), NominalZd.GetArray(), n);
    fSplineAz = new TSpline3("azimuth",
			     ReportTime.GetArray(), NominalAz.GetArray(), n);
    fSplineRa = new TSpline3("RA",
			     ReportTime.GetArray(), Ra.GetArray(), n);
    fSplineDec = new TSpline3("DEC",
			     ReportTime.GetArray(), Dec.GetArray(), n);


  //    TCanvas* c = new TCanvas();
//      c->Divide(2,2);
//      c->cd(1);
//      fSplineZd->Draw();
//      c->cd(2);
//      fSplineAz->Draw();
//      c->cd(3);
//      fSplineRa->Draw(); 
//      c->cd(4);
//      fSplineDec->Draw();
}


// --------------------------------------------------------------------------
//
// default constructor
//
MExtrapolatePointingPos::MExtrapolatePointingPos(const TString filename, const char *name, const char *title)
    
{
    fName  = name  ? name  : "MExtrapolatePointingPos";
    fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";

   
    fFilename = filename;
}

MExtrapolatePointingPos::~MExtrapolatePointingPos()
{
    delete fSplineZd;
    delete fSplineAz;
    delete fSplineRa;    
    delete fSplineDec;
}


// --------------------------------------------------------------------------
//
//  Input:
//  - MTime
// 
//  Output:
//  - MPointingPos
//
Int_t MExtrapolatePointingPos::PreProcess( MParList *pList )
{
    fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
    if (!fRunHeader)
    {
        *fLog << err << "MRunHeader not found... aborting." << endl;
        return kFALSE;
    }
    
    fEvtTime = (MTime*)pList->FindObject("MTime");
    if (!fEvtTime)
    {
        *fLog << err << "MTime not found... aborting." << endl;
        return kFALSE;
    }

    fPointingPos = (MPointingPos*)pList->FindCreateObj("MPointingPos");
    if (!fPointingPos)
        return kFALSE;


    ReadDriveReport(fFilename);


    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Fill the MPedestalCam container with the signal mean and rms for the event.
// Store the measured signal in arrays fSumx and fSumx2 so that we can 
// calculate the overall mean and rms in the PostProcess()
//
Int_t MExtrapolatePointingPos::Process()
{

    const MTime* StartRunTime  = &fRunHeader->GetRunStart();
    Int_t run = fRunHeader->GetRunNumber();

    Long_t time;

     //   if(run < 20000)
//    	time = fEvtTime->GetTime();
//        else
	time = StartRunTime->GetTime();


 
    
    Double_t zd = fSplineZd->Eval( time );
    Double_t az = fSplineAz->Eval( time );
    Double_t ra = fSplineRa->Eval( time );
    Double_t dec = fSplineDec->Eval( time );

    fPointingPos->SetLocalPosition( zd, az );
    fPointingPos->SetSkyPosition( ra*TMath::DegToRad()/15, dec*TMath::DegToRad());

    //  *fLog << " PointingPos: " << " time = " << time << " " << *fEvtTime << " (zd, az, ra, dec) = (" << zd << ", "  << az << ", " << ra << ", "  << dec << ")" << endl; 


  return kTRUE;
}


