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

/* ======================================================================== *\
!
! *
! * 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 "MPointingDev.h"
#include "MRawRunHeader.h"
#include "MReportStarguider.h"

ClassImp(MPointingDevCalc);

using namespace std;

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

    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;

    fSkip.Reset();

    return fDeviation ? kTRUE : kFALSE;
}

Int_t MPointingDevCalc::ProcessStarguiderReport()
{
    Double_t devzd = fReport->GetDevZd(); // [deg]
    Double_t devaz = fReport->GetDevAz(); // [deg]
    if (devzd==0 && devaz==0)
    {
        fDeviation->SetDevZdAz(0, 0);
        fSkip[1]++;
        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)
            {
                fSkip[2]++;
                return kTRUE;
            }
        }

        if (fReport->GetNumIdentifiedStars()<fNumMinStars)
        {
            fSkip[3]++;
            return kTRUE;
        }
    }

    // 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)
    {
        fSkip[4]++;
        return kTRUE;
    }

    fDeviation->SetDevZdAz(devzd, devaz);

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

    fSkip[0]++;

    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);
        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], Form("NSB out of %.1f sigma range", fNsbLevel));
    PrintSkipped(fSkip[3], Form("Number of identified stars < %d", fNumMinStars));
    PrintSkipped(fSkip[4], Form("Absolute deviation > %.1farcmin", fMaxAbsDev));
    *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
//
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;
    }

    return rc;
}