source: trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc@ 8601

/* ======================================================================== *\
!
! *
! * 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, 07/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MPointingDevCalc
//
// Calculates the pointing deviation from the starguider information.
//
// There are some quality parameters to steer which are the valid
// starguider data points:
//  * MPointingDevCalc.NumMinStars: 8
//    Minimum number of identified stars required to accep the data
//  * MPointingDevCalc.NsbLevel:    3.0
//    Minimum deviation (in rms) from the the mean allowed for the measured
//    NSB (noise in the ccd camera)
//  * MPointingDevCalc.NsbMin:      30
//    - minimum NSB to be used in mean/rms calculation
//  * MPointingDevCalc.NsbMax:      60
//    - maximum NSB to be used in mean/rms calculation
//  * MPointingDevCalc.MaxAbsDev:   15
//    - Maximum absolute deviation which is consideres as valid (arcmin)
//
// Starguider data which doens't fullfill this requirements is ignored.
// If the measures NSB==0 (file too old, starguider didn't write down
// these values) the checks based on NSB and NumMinStar are skipped.
//
// The calculation of NSB mean and rms is reset for each file (ReInit)
//
//
// Input Container:
//   MReportStarguider
//
// Output Container:
//   MPointingDev
//
/////////////////////////////////////////////////////////////////////////////
#include "MPointingDevCalc.h"

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

#include "MParList.h"

#include "MAstro.h"
#include "MPointing.h"
#include "MPointingDev.h"
#include "MRawRunHeader.h"
#include "MReportStarguider.h"

ClassImp(MPointingDevCalc);

using namespace std;

const TString MPointingDevCalc::fgFileName="resources/starguider.txt";

// --------------------------------------------------------------------------
//
// Delete fPointing and set pointing to NULL
//
void MPointingDevCalc::Clear(Option_t *o)
{
    if (fPointing)
        delete fPointing;

    fPointing = NULL;
}

// --------------------------------------------------------------------------
//
// Clear the pointing model. If run-number >= 87751 read the new
// pointing model with fFileName
//
Bool_t MPointingDevCalc::ReadPointingModel(const MRawRunHeader &run)
{
    if (run.GetRunNumber()<87751)
    {
        Clear();
        return kTRUE;
    }

    if (!fPointing)
        fPointing = new MPointing;

    if (fFileName==fPointing->GetName())
    {
        *fLog << inf << fFileName << " already loaded." << endl;
        return kTRUE;
    }

    return fPointing->Load(fFileName);
}

// --------------------------------------------------------------------------
//
// Check the file/run type from the run-header and if it is a data file
// load starguider calibration.
//
Bool_t MPointingDevCalc::ReInit(MParList *plist)
{
    MRawRunHeader *run = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
    if (!run)
    {
        *fLog << err << "MRawRunHeader not found... aborting." << endl;
        return kFALSE;
    }

    fNsbSum   =  0;
    fNsbSq    =  0;
    fNsbCount =  0;

    fRunType = run->GetRunType();

    switch (fRunType)
    {
    case MRawRunHeader::kRTData:
        if (!fReport)
            *fLog << warn << "MReportStarguider not found... skipped." << endl;
        return ReadPointingModel(*run);

    case MRawRunHeader::kRTMonteCarlo:
        return kTRUE;

    case MRawRunHeader::kRTPedestal:
        *fLog << err << "Cannot work in a pedestal Run!... aborting." << endl;
        return kFALSE;

    case MRawRunHeader::kRTCalibration:
        *fLog << err << "Cannot work in a calibration Run!... aborting." << endl;
        return kFALSE;

    default:
        *fLog << err << "Run Type " << fRunType << " unknown!... aborting." << endl;
        return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Search for 'MPointingPos'. Create if not found.
//
Int_t MPointingDevCalc::PreProcess(MParList *plist)
{
    fDeviation = (MPointingDev*)plist->FindCreateObj("MPointingDev");
    fReport    = (MReportStarguider*)plist->FindObject("MReportStarguider");

    // We use kRTNone here as a placeholder for data runs.
    fRunType  = MRawRunHeader::kRTNone;
    fLastMjd  = -1;

    fSkip.Reset();

    return fDeviation ? kTRUE : kFALSE;
}

// --------------------------------------------------------------------------
//
// Increase fSkip[i] by one. If the data in fDeviation is outdated (older
// than fMaxAge) and the current report should be skipped reset DevZdAz and
// DevXY and fSkip[6] is increased by one.
//
void MPointingDevCalc::Skip(Int_t i)
{
    fSkip[i]++;

    const Double_t diff = (fReport->GetMjd()-fLastMjd)*1440; // [min] 1440=24*60
    if (diff<fMaxAge && fLastMjd>0)
        return;

    fDeviation->SetDevZdAz(0, 0);
    fDeviation->SetDevXY(0, 0);
    fSkip[6]++;
}

// --------------------------------------------------------------------------
//
// Do a full starguider calibration using a pointing model for the starguider.
//
void MPointingDevCalc::DoCalibration(Double_t devzd, Double_t devaz) const
{
    if (!fPointing)
    {
        // Do a simple starguider calibration using a simple offset in x and y
        fDeviation->SetDevZdAz(devzd, devaz);

        // Linear starguider calibration taken from April/May data
        // For calibration add MDriveReport::GetErrorZd/Az !
        fDeviation->SetDevXY(fDx, fDy); // 1arcmin ~ 5mm
        //devzd -= 2.686/60;   
        //devaz -= 2.840/60;

        return;
    }

    // def?: 20.105, 773
    // 0/0 : 20.119, 763
    // -/- : 19.417  726
    // +mis: 19.80,  756

    // Get the nominal position the star is at the sky
    // Unit: deg
    ZdAz nom(fReport->GetNominalZd(), fReport->GetNominalAz());
    nom *= TMath::DegToRad();

    // Get the mispointing measured by the telescope. It is
    // calculate as follows:
    //
    // Position at which the starguider camera is pointing in real:
    //       pointing position = nominal position - mis
    ZdAz mis(devzd, devaz);
    mis *= TMath::DegToRad();

    // The pointing mode is the conversion from the real pointing
    // position of the telescope into the pointing position measured
    // by the starguider.
    //
    //  --> To get the real poiting position of the telescope we have
    //      to convert the measured position back;
    //

    // Position at which the starguider camera is pointing in real:
    //       pointing position = nominal position - dev
    //
    // The position measured as the starguider's pointing position
    ZdAz pos(nom);        // cpos = sao - dev
    pos -= mis;

    // Now we convert the starguider's pointing position into the
    // telescope pointing position (the pointing model converts
    // the telescope pointing position into the starguider pointing
    // position)
    ZdAz point = fPointing->CorrectBack(pos);

    // MSrcPosCalc uses the following condition to calculate the
    // source position in the camera:
    //    real pointing pos = nominal pointing pos - dev
    //
    // Therefor we calculate dev as follows:
    ZdAz dev(nom);
    dev -= point;
    dev *= TMath::RadToDeg();

    // Set Measured mispointing and additional offset
    fDeviation->SetDevZdAz(dev.Zd(), dev.Az());
    fDeviation->SetDevXY(0, 0);
}

Int_t MPointingDevCalc::ProcessStarguiderReport()
{
    /************* CHECK STATUS!!! ******************/

    Double_t devzd = fReport->GetDevZd(); // [arcmin]
    Double_t devaz = fReport->GetDevAz(); // [arcmin]
    if (devzd==0 && devaz==0)
    {
        Skip(1);
        return kTRUE;
    }

    if (!fReport->IsMonitoring())
    {
        Skip(2);
        return kTRUE;
    }

    devzd /= 60; // Convert from arcmin to deg
    devaz /= 60; // Convert from arcmin to deg

    const Double_t nsb = fReport->GetSkyBrightness();
    if (nsb>0)
    {
        if (nsb>fNsbMin && nsb<fNsbMax)
        {
            fNsbSum += nsb;
            fNsbSq  += nsb*nsb;
            fNsbCount++;
        }

        if (fNsbCount>0)
        {
            const Double_t sum = fNsbSum/fNsbCount;
            const Double_t sq  = fNsbSq /fNsbCount;

            const Double_t rms = fNsbLevel*TMath::Sqrt(sq - sum*sum);

            if (nsb<sum-rms || nsb>sum+rms)
            {
                Skip(3);
                return kTRUE;
            }
        }

        if (fReport->GetNumIdentifiedStars()<fNumMinStars)
        {
            Skip(4);
            return kTRUE;
        }
    }

    // >= 87751 (31.3.06 12:00)

    // Calculate absolute deviation
    const Double_t dev = MAstro::GetDevAbs(fReport->GetNominalZd(), devzd, devaz);

    // Sanity check... larger deviation are strange and ignored
    if (dev*60>fMaxAbsDev)
    {
        Skip(5);
        return kTRUE;
    }

    DoCalibration(devzd, devaz);

    fSkip[0]++;
    fLastMjd = fReport->GetMjd();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  See class description.
//
Int_t MPointingDevCalc::Process()
{
    switch (fRunType)
    {
    case MRawRunHeader::kRTNone:
    case MRawRunHeader::kRTData:
        return fReport ? ProcessStarguiderReport() : kTRUE;

    case MRawRunHeader::kRTMonteCarlo:
        fSkip[0]++;
        fDeviation->SetDevZdAz(0, 0);
        fDeviation->SetDevXY(0, 0);
        return kTRUE;
    }
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Print execution statistics
//
Int_t MPointingDevCalc::PostProcess()
{
    if (GetNumExecutions()==0)
        return kTRUE;

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;
    PrintSkipped(fSkip[1], "Starguider deviation not set, is exactly 0/0");
    PrintSkipped(fSkip[2], "Starguider was not monitoring (eg. LEDs off)");
    PrintSkipped(fSkip[3], Form("NSB out of %.1f sigma range", fNsbLevel));
    PrintSkipped(fSkip[4], Form("Number of identified stars < %d", fNumMinStars));
    PrintSkipped(fSkip[5], Form("Absolute deviation > %.1farcmin", fMaxAbsDev));
    PrintSkipped(fSkip[6], Form("Events set to 0 because older than %.1fmin", fMaxAge));
    *fLog << " " << (int)fSkip[0] << " (" << Form("%5.1f", 100.*fSkip[0]/GetNumExecutions()) << "%) Evts survived calculation!" << endl;
    *fLog << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// MPointingDevCalc.NumMinStars: 8
// MPointingDevCalc.NsbLevel:    3.0
// MPointingDevCalc.NsbMin:      30
// MPointingDevCalc.NsbMax:      60
// MPointingDevCalc.MaxAbsDev:   15
// MPointingDevCalc.MaxAge:      1.0
//
Int_t MPointingDevCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "NumMinStars", print))
    {
        SetNumMinStars(GetEnvValue(env, prefix, "NumMinStars", (Int_t)fNumMinStars));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "NsbLevel", print))
    {
        SetNsbLevel(GetEnvValue(env, prefix, "NsbLevel", fNsbLevel));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "NsbMin", print))
    {
        SetNsbMin(GetEnvValue(env, prefix, "NsbMin", fNsbMin));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "NsbMax", print))
    {
        SetNsbMax(GetEnvValue(env, prefix, "NsbMax", fNsbMax));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "MaxAbsDev", print))
    {
        SetMaxAbsDev(GetEnvValue(env, prefix, "MaxAbsDev", fMaxAbsDev));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "Dx", print))
    {
        fDx = GetEnvValue(env, prefix, "Dx", fDx);
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "Dy", print))
    {
        fDy = GetEnvValue(env, prefix, "Dy", fDy);
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "MaxAge", print))
    {
        fMaxAge = GetEnvValue(env, prefix, "MaxAge", fMaxAge);
        rc = kTRUE;
    }

    return rc;
}