source: trunk/MagicSoft/Mars/mtemp/mucm/classes/MExtrapolatePointingPos.cc@ 4691

/* ======================================================================== *\
!
! *
! * 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;


// ---------------------------------------------------------------------------
//
// Read the drive report file for the whole night, a build from it the splines
//
Bool_t MExtrapolatePointingPos::ReadDriveReport(const TString filename)
{

    *fLog << endl << "["<< GetName() << "]: Loading report file \"" << filename << "\" into TSpline..." << endl;


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


    //
    // TaskList
    //
    MReadReports read;
    read.AddTree("Drive");
    read.AddFile(filename);     // after the reading of the trees!!!
    read.AddToBranchList("MReportDrive.*");
    
    tlist.AddToList(&read);
 
   
    //
    // EventLoop
    //
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);

    if (!evtloop.PreProcess())
        return kFALSE;


    TArrayD ReportTime(10000);
    TArrayD CurrentZd(10000);
    TArrayD CurrentAz(10000); 
    TArrayD NominalZd(10000);
    TArrayD NominalAz(10000);
    TArrayD Ra(10000);
    TArrayD Dec(10000);

    
    Int_t n=0;
    while (tlist.Process())
    { 
        MReportDrive* report = (MReportDrive*)plist.FindObject("MReportDrive");
        MTime* reporttime = (MTime*)plist.FindObject("MTimeDrive");
        
        if(n==0)
            fFirstDriveTime = *reporttime;
        else 
            fLastDriveTime = *reporttime;

        //
        // Sometimes there are two reports with the same time
        //
        if (reporttime->GetTime() == ReportTime[n-1])
        { 
            cout << warn <<"["<< GetName() << "]: Warning: this report has the same time that the previous one...skipping it " << 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();

        n++;
    }

    tlist.PrintStatistics();

    *fLog << "["<< GetName() << "]: loaded " << n << " ReportDrive from " 
          << fFirstDriveTime << " to " << fLastDriveTime << endl << endl;


    //
    // 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);
     
    // for(int i=0;i<ReportTime.GetSize();i++)
//      *fLog << i            << " " << ReportTime[i] << " " 
//            << CurrentZd[i] << " " << CurrentAz[i]  << " " 
//            << NominalZd[i] << " " << NominalAz[i]  << " " 
//            << Ra[i]        << " " << Dec[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();
//     c->Modified();
//     c->Update();
    
    return kTRUE;
}


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

    // Init
    fSplineZd  = NULL;
    fSplineAz  = NULL;
    fSplineRa  = NULL;
    fSplineDec = NULL;
}

MExtrapolatePointingPos::~MExtrapolatePointingPos()
{
    if(fSplineZd)
        delete fSplineZd;
    if(fSplineAz)
        delete fSplineAz;
    if(fSplineRa)
        delete fSplineRa;    
    if(fSplineDec)
        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;


    if( !ReadDriveReport(fFilename) )
        return kFALSE;


    return kTRUE;
}


// --------------------------------------------------------------------------
//
//  Get the run start time, and get the pointing position for that time
//
Int_t MExtrapolatePointingPos::Process()
{

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


    //
    // Check that we have drive report for this time
    //
    if( *StartRunTime<fFirstDriveTime || *StartRunTime>fLastDriveTime)
    {
        *fLog << err << dbginf << GetName() << ": Run time " << *StartRunTime
              << " outside range of drive reports  (" << fFirstDriveTime 
              << ", " << fLastDriveTime << ")" << endl;
        return kFALSE;
    }

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


    const Double_t zd = fSplineZd->Eval( time );
    const Double_t az = fSplineAz->Eval( time );
    const Double_t ra = fSplineRa->Eval( time );
    const 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;
}