source: trunk/MagicSoft/Mars/mastro/MAstroCamera.cc@ 3537

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

/////////////////////////////////////////////////////////////////////////////
//
// MAstroCamera
//
// A tools displaying stars from a catalog in the camera display.
//
// For a usage example see macros/starfield.C
//
// PRELIMINARY!!
//
/////////////////////////////////////////////////////////////////////////////
#include "MAstroCamera.h"

#include <KeySymbols.h>

#include <TH2.h>
#include <TMarker.h>
#include <TVirtualPad.h>

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

#include "MGeomCam.h"
#include "MGeomMirror.h"

#include "MTime.h"
#include "MAstroSky2Local.h"
#include "../mhist/MHCamera.h"
#include "MObservatory.h"

ClassImp(MAstroCamera);

using namespace std;

MAstroCamera::MAstroCamera() : fGeom(0), fMirrors(0)
{
    fMirror0 = new MGeomMirror;
    fMirror0->SetMirrorContent(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1);
}

MAstroCamera::~MAstroCamera()
{
    if (fGeom)
        delete fGeom;
    if (fMirrors)
        delete fMirrors;

    delete fMirror0;
}

void MAstroCamera::SetMirrors(TClonesArray *arr)
{
    if (!arr || arr->GetClass()!=MGeomMirror::Class())
        return;

    const Int_t n = arr->GetSize();

    if (!fMirrors)
        fMirrors = new TClonesArray(MGeomMirror::Class(), n);

    fMirrors->ExpandCreate(n);

    for (int i=0; i<n; i++)
        memcpy((*fMirrors)[i], (*arr)[i], sizeof(MGeomMirror));

}

void MAstroCamera::SetGeom(const MGeomCam &cam)
{
    if (fGeom)
        delete fGeom;

    fGeom = (MGeomCam*)cam.Clone();
}

Int_t MAstroCamera::ConvertToPad(const TVector3 &w, TVector2 &v)
{
    const TVector3 spot = fMirror0->GetReflection(w, fGeom->GetCameraDist())*1000;

    /*
     --- Use this to plot the 'mean grid' instead of the grid of a
         theoretical central mirror ---

        TVector3 spot;
        const Int_t num = fConfig->GetNumMirror();
        for (int i=0; i<num; i++)
        spot += fConfig->GetMirror(i).GetReflection(w, fGeom->GetCameraDist())*1000;
        spot *= 1./num;
        */


    v.Set(spot(0), spot(1));

    const Float_t max = fGeom->GetMaxRadius()*0.70;
    return v.X()>-max && v.Y()>-max && v.X()<max && v.Y()<max;
}

void MAstroCamera::DrawNet(const TRotation &trans)
{
    const TRotation rot(MAstroSky2Local(*fTime, *fObservatory));

    TVector2 radec(fRaDec.Phi(), fRaDec.Theta());
    MAstroCatalog::DrawNet(radec, trans*rot, 2);

    const TVector3 zdaz0 = MAstroSky2Local(*fTime, *fObservatory)*fRaDec;
    TVector2 zdaz(zdaz0.Phi(), zdaz0.Theta());
    MAstroCatalog::DrawNet(zdaz, trans, 1);
}

TObject *FindObjectInPad(const char *name, TVirtualPad *pad)
{
    if (!pad)
        pad = gPad;

    if (!pad)
        return NULL;

    TObject *o;

    TIter Next(pad->GetListOfPrimitives());
    while ((o=Next()))
    {
        if (o->InheritsFrom(gROOT->GetClass(name)))
            return o;

        if (o->InheritsFrom("TPad"))
            if ((o = FindObjectInPad(name, (TVirtualPad*)o)))
                return o;
    }
    return NULL;
}

void MAstroCamera::AddPrimitives(Option_t *o)
{
    if (!fTime || !fObservatory || !fMirrors)
    {
        cout << "Missing data..." << endl;
        return;
    }

    TString opt(o);
    if (opt.IsNull())
        opt = "*.";

    const Bool_t hashist = opt.Contains("h", TString::kIgnoreCase);
    const Bool_t hasmean = opt.Contains("*", TString::kIgnoreCase);
    const Bool_t hasdot  = opt.Contains(".", TString::kIgnoreCase);
    const Bool_t usecam  = opt.Contains("c", TString::kIgnoreCase);

    MAstroSky2Local rot(*fTime, *fObservatory);

    const Float_t rho = rot.RotationAngle(fRaDec.Phi(), TMath::Pi()/2-fRaDec.Theta());

    TString str = fTime->GetSqlDateTime();
    str += Form("  (\\alpha=%.1fh \\delta=%.1f\\circ)  \\rho=%.1f\\circ",
                fRaDec.Phi()/TMath::Pi()*12, 90-fRaDec.Theta()*TMath::RadToDeg(),
                rho *TMath::RadToDeg());

    // Get camera
    MHCamera *camera=(MHCamera*)FindObjectInPad("MHCamera", gPad);
    if (camera)
    {
        if (!camera->GetGeometry() || camera->GetGeometry()->IsA()!=fGeom->IsA())
            camera->SetGeometry(*fGeom);
    }
    else
    {
        camera = new MHCamera(*fGeom);
        camera->SetName("MHCamera");
        camera->SetStats(0);
        camera->SetInvDeepBlueSeaPalette();
        camera->SetBit(kCanDelete);
        camera->Draw();
    }

    camera->SetTitle(str);

    gPad->cd(1);

    if (!usecam)
    {
        if (camera->GetEntries()==0)
            camera->SetBit(MHCamera::kNoLegend);
    }
    else
    {
        camera->Reset();
        camera->SetYTitle("arb.cur");
    }

    TH2 *h=0;
    if (hashist)
    {
        TH2F hist("","", 90, -650, 650, 90, -650, 650);
        hist.SetMinimum(0);
        h = (TH2*)hist.DrawCopy("samecont1");
    }

    TVector3 zdaz0 = fRaDec;
    zdaz0 *= rot;

    cout << zdaz0.Phi()*TMath::RadToDeg() << " " << zdaz0.Theta()*TMath::RadToDeg() << endl;

    TVector3 test = zdaz0;
    test *= rot.Inverse();

    cout << test.Phi()*TMath::RadToDeg()/15 << " " << test.Theta()*TMath::RadToDeg() << endl;

    TRotation rot2;
    rot2.RotateZ(-zdaz0.Phi());
    rot2.RotateY(-zdaz0.Theta());
    rot2.RotateZ(-TMath::Pi()/2); // align coordinate system

    DrawNet(rot2);

    MVector3 *radec;
    TIter Next(&fList);

    while ((radec=(MVector3*)Next()))
    {
        const Double_t mag = radec->Magnitude();

        TVector3 star(*radec);

        // Calculate local coordinates
        star *= rot;
        // Rotate Star into telescope system
        star *= rot2;

        TVector3 mean;

        Int_t num = 0;

        MGeomMirror *mirror = 0;
        TIter NextM(fMirrors);
        while ((mirror=(MGeomMirror*)NextM()))
        {
            const TVector3 spot = mirror->GetReflection(star, fGeom->GetCameraDist())*1000;

            // calculate mean of all 'stars' hitting the camera plane
            // by taking the sum of the intersection points between
            // the light vector and the camera plane
            mean += spot;

            if (hasdot)
            {
                TMarker *m=new TMarker(spot(0), spot(1), 1);
                m->SetBit(kCannotPick);
                m->SetBit(kCanDelete);
                m->SetMarkerColor(kMagenta);
                m->SetMarkerStyle(kDot);
                AddMap(m);
            }
            if (h)
                h->Fill(spot(0), spot(1), pow(10, -mag/2.5));

            if (usecam)
                camera->Fill(spot(0), spot(1), pow(10, -mag/2.5));

            num++;
        }

        // transform meters into millimeters (camera display works with mm)
        mean *= 1./num;

        DrawStar(mean(0), mean(1), *radec, !hasmean, Form("x=%.1fmm y=%.1fmm", mean(0), mean(1)));
    }
}

// ------------------------------------------------------------------------
//
// Execute a gui event on the camera
//
void MAstroCamera::ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2)
{
    if (event==kKeyPress && fTime)
        switch (mp2)
        {
        case kKey_Plus:
            fTime->SetMjd(fTime->GetMjd()+0.25/24);
            SetBit(kHasChanged);
            gPad->Modified();
            gPad->Update();
            return;

        case kKey_Minus:
            fTime->SetMjd(fTime->GetMjd()-0.25/24);
            SetBit(kHasChanged);
            gPad->Modified();
            gPad->Update();
            return;
        }

    MAstroCatalog::ExecuteEvent(event, mp1, mp2);
}