source: trunk/MagicSoft/Mars/mimage/MHillasSrc.cc@ 4256

/* ======================================================================== *\
!
! *
! * 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    12/2000 <mailto:tbretz@uni-sw.gwdg.de>
!   Author(s): Harald Kornmayer 1/2001
!   Author(s): Rudolf Bock     10/2001 <mailto:Rudolf.Bock@cern.ch>
!   Author(s): Wolfgang Wittek 06/2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHillasSrc
//
// Storage Container for image parameters
//
//    source-dependent image parameters
//
//
// Version 1:
// ----------
//  fAlpha          angle between major axis and line source-to-center
//  fDist           distance from source to center of ellipse
//
//
// Version 2:
// ----------
//  fHeadTail       added
//
//
// Version 3:
// ----------
//  fCosDeltaAlpha  cosine of angle between d and a, where
//                   - d is the vector from the source position to the
//                     center of the ellipse
//                   - a is a vector along the main axis of the ellipse,
//                     defined with positive x-component
//
//
// Version 4:
// ----------
//
// fHeadTail        removed
//
/////////////////////////////////////////////////////////////////////////////
#include "MHillasSrc.h"

#include <fstream>
#include <TArrayF.h>

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

#include "MGeomCam.h"
#include "MSrcPosCam.h"

ClassImp(MHillasSrc);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MHillasSrc::MHillasSrc(const char *name, const char *title)
{
    fName  = name  ? name  : "MHillasSrc";
    fTitle = title ? title : "Parameters depending in source position";
}

void MHillasSrc::Reset()
{
    fDist          = -1;
    fAlpha         =  0;
    fCosDeltaAlpha =  0;
}

// --------------------------------------------------------------------------
//
//  Calculation of source-dependent parameters
//  In case you don't call Calc from within an eventloop make sure, that
//  you call the Reset member function before.
//
Bool_t MHillasSrc::Calc(const MHillas *hillas)
{
    const Double_t mx = hillas->GetMeanX();     // [mm]
    const Double_t my = hillas->GetMeanY();     // [mm]

    const Double_t sx = mx - fSrcPos->GetX();   // [mm]
    const Double_t sy = my - fSrcPos->GetY();   // [mm]

    const Double_t sd = hillas->GetSinDelta();  // [1]
    const Double_t cd = hillas->GetCosDelta();  // [1]

    //
    // Distance from source position to center of ellipse.
    // If the distance is 0 distance, Alpha is not specified.
    // The calculation has failed and returnes kFALSE.
    //
    const Double_t dist = sqrt(sx*sx + sy*sy);  // [mm]
    if (dist==0)
        return kFALSE;

    //
    // Calculate Alpha and Cosda = cos(d,a)
    // The sign of Cosda will be used for quantities containing 
    // a head-tail information
    //
    // *OLD* const Double_t arg = (sy-tand*sx) / (dist*sqrt(tand*tand+1));
    // *OLD* fAlpha = asin(arg)*kRad2Deg;
    //
    const Double_t arg1 = cd*sy-sd*sx;          // [mm]
    const Double_t arg2 = cd*sx+sd*sy;          // [mm]

    //
    // Due to numerical uncertanties in the calculation of the
    // square root (dist) and arg1 it can happen (in less than 1e-5 cases)
    // that the absolute value of arg exceeds 1. Because this uncertainty
    // results in an Delta Alpha which is still less than 1e-3 we don't care
    // about this uncertainty in general and simply set values which exceed
    // to 1 saving its sign.
    //
    const Double_t arg = arg1/dist;
    fAlpha = TMath::Abs(arg)>1 ? TMath::Sign(90., arg) : asin(arg)*kRad2Deg;  // [deg]

    fCosDeltaAlpha = arg2/dist;                 // [1]
    fDist          = dist;                      // [mm]

    SetReadyToSave();

    return kTRUE;
} 

// --------------------------------------------------------------------------
//
// Print contents of MHillasSrc to *fLog
//
void MHillasSrc::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "Source dependant Image Parameters (" << GetName() << ")" << endl;
    *fLog << " - Dist           [mm]  = " << fDist << endl;
    *fLog << " - Alpha          [deg] = " << fAlpha << endl;
    *fLog << " - CosDeltaAlpha        = " << fCosDeltaAlpha << endl;
}

// --------------------------------------------------------------------------
//
// Print contents of MHillasSrc to *fLog depending on the geometry in
// units of deg.
//
void MHillasSrc::Print(const MGeomCam &geom) const
{
    *fLog << all;
    *fLog << "Source dependant Image Parameters (" << GetName() << ")" << endl;
    *fLog << " - Dist           [deg] = " << fDist*geom.GetConvMm2Deg() << endl;
    *fLog << " - Alpha          [deg] = " << fAlpha << endl;
    *fLog << " - CosDeltaAlpha        = " << fCosDeltaAlpha << endl;
}

// --------------------------------------------------------------------------
//
// This function is ment for special usage, please never try to set
// values via this function
//
void MHillasSrc::Set(const TArrayF &arr)
{
    if (arr.GetSize() != 3)
        return;

    fAlpha = arr.At(0);         // [deg]  angle of major axis with vector to src
    fDist  = arr.At(1);         // [mm]   distance between src and center of ellipse
    fCosDeltaAlpha = arr.At(2); // [1]    cosine of angle between d and a
}