source: trunk/Cosy/tpoint/plot.C@ 10075

#include <iostream>
#include <iomanip>
#include <fstream>

#include <TGraph.h>
#include <TCanvas.h>
#include <TSystem.h>

#include "mars/MTime.h"
#include "mars/MMath.h"
#include "mars/MAstro.h"
#include "mars/MDirIter.h"
#include "mars/MStatusDisplay.h"

using namespace std;

int ReadFile(const char *fname, TGraph &g1, TGraph &g2, TGraph &g0)
{
    ifstream fin(fname);

    cout << "Reading " << setw(23) << gSystem->BaseName(fname) << "..." << flush;

    MTime t;
    while (1)
    {
        TString str;
        str.ReadLine(fin);
        if (!fin)
            break;

        if (str.Contains("#"))
            continue;

        Double_t alt, az, dalt, daz, mjd;
        if (sscanf(str.Data(), "%lf %lf %*f %*f %*f %*f %lf %lf %lf",
                   &az, &alt, &dalt, &daz, &mjd)!=5)
            continue;

        if (dalt==0)
            continue;

        Double_t res = MAstro::GetDevAbs(90-(alt-dalt), -dalt, -daz)*60;

        if (res>12)
            continue;

        t.SetMjd(mjd);

        g0.SetPoint(g0.GetN(), t.GetAxisTime(), g0.GetN());
        g1.SetPoint(g1.GetN(), t.GetAxisTime(), res);
        g2.SetPoint(g2.GetN(), g2.GetN(),       res);

        if (g0.GetN()<2 || g0.GetX()[g0.GetN()-2]<g0.GetX()[g0.GetN()-1])
            continue;

        cout << "WARNING-TIME-SKEW! " << mjd << " " << flush;
    }

    cout << g1.GetN() << endl;

    return g1.GetN();
}

void DrawGrid(TH1 &h)
{
    TGraph l;
    l.SetLineStyle(kDashed);
    l.SetLineWidth(2);
    l.SetLineColor(kBlue);

    Int_t style[] = { kSolid, 9, 7 };

    Double_t x[2] = { h.GetXaxis()->GetXmin(), h.GetXaxis()->GetXmax() };
    Double_t y[2];

    for (int i=0; i<3; i++)
    {
        y[0] = y[1] = i+1;
        l.SetLineStyle(style[i]);
        l.DrawGraph(2, x, y, "L");
    }

    if (x[0]<1e6)
        return;

    Double_t X[2];
    Double_t Y[2] = { h.GetMinimum(), h.GetMaximum() };

    l.SetLineWidth(1);

    for (int i=2004; i<2020; i++)
    {
        for (int j=0; j<12; j++)
        {
            X[0] = X[1] = MTime(i, 1+j, 1).GetAxisTime();

            if (X[0]<x[0] || X[0]>x[1])
                continue;

            l.SetLineStyle(j==0 ? kDashed : kDotted);
            l.DrawGraph(2, X, Y, "L");
        }
    }
}

void DrawTimes(TH1 &h)
{
    TGraph l;
    l.SetLineColor(kBlue);

    Double_t x[2];
    Double_t y[2] = { h.GetMinimum(), h.GetMaximum() };

    for (int i=0; dates[i+1]>0; i++)
    {
        x[0] = x[1] = dates[i];
        l.DrawGraph(2, x, y, "L");
    }

    DrawGrid(h);
}

Int_t GetBin(TGraph &times, Int_t i)
{
    return TMath::BinarySearch(times.GetN(), times.GetX(), dates[i]);;
}

void DrawIndices(TGraph &times, TH1 &h)
{
    TGraph l;
    l.SetLineStyle(kDashed);
    l.SetLineColor(kBlue);

    Double_t x[2];
    Double_t y[2] = { h.GetMinimum(), h.GetMaximum() };

    for (int i=0; dates[i+1]>0; i++)
    {
        Int_t bin = GetBin(times, i);
        if (bin<0)
            continue;

        x[0] = x[1] = bin+0.5;
        l.DrawGraph(2, x, y, "L");
    }

    DrawGrid(h);
}

void plot(MDirIter &Next)
{
    Next.Sort();

    TGraph g1;
    TGraph g2;
    TGraph g0;

    while (1)
    {
        TString name=Next();
        if (name.IsNull())
            break;

        ReadFile(name, g1, g2, g0);
    }

    MStatusDisplay *d = new MStatusDisplay;

    gROOT->SetSelectedPad(0);
    TCanvas &c0 = d->AddTab("TimeLine");
    c0.SetGrid();

    g0.SetTitle("TimeLine");
    g0.GetXaxis()->SetTimeDisplay(1);
    g0.SetMarkerStyle(kFullDotMedium);
    g0.SetLineWidth(1);
    g0.SetLineColor(kGreen);
    g0.DrawClone("ALP");

    DrawTimes(*g0.GetHistogram());


    gROOT->SetSelectedPad(0);
    TCanvas &c1 = d->AddTab("ResDate");
    //c1.SetGrid();

    g1.SetTitle("Measured Residual vs. Date");
    g1.GetXaxis()->SetTimeDisplay(1);
//    g1.SetMaximum(10);
    g1.SetMarkerStyle(kFullDotMedium);
    g1.DrawClone("AP");

    DrawTimes(*g1.GetHistogram());


    gROOT->SetSelectedPad(0);
    TCanvas &c2 = d->AddTab("ResIdx");
    c2.SetGrid();

    g2.SetTitle("Measured Residual vs. Index");
//    g2.SetMaximum(0.1);
    g2.SetMarkerStyle(kFullDotMedium);
    g2.DrawClone("AP");

    DrawIndices(g0, *g2.GetHistogram());

    Double_t prob[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
    Double_t quant[4];

    int n;
    for (n=0; dates[n]>0; n++);

    TH1F h0("Q98", "", n-1, dates);
    TH1F h1("Q95", "", n-1, dates);
    TH1F h2("Q68", "", n-1, dates);
    TH1F h3("Q50", "", n-1, dates);

    Int_t first = 0;

    cout << "n=" << n << endl;

    for (int i=0; dates[i+1]>0; i++)
    {
        Int_t bin0 = GetBin(g0, i);
        Int_t bin1 = GetBin(g0, i+1);

        cout << i << " " << bin0 << " " << bin1 << endl;
        if (bin0<0)
            continue;

        if (bin0==bin1)
            continue;

        TMath::Quantiles(bin1-bin0, 4, g2.GetY()+bin0, quant, prob, kFALSE, NULL, 7);

        h0.SetBinContent(i+1, quant[0]);
        h1.SetBinContent(i+1, quant[1]);
        h2.SetBinContent(i+1, quant[2]);
        h3.SetBinContent(i+1, quant[3]);

        if (first==0)
            first = i;
    }


    gROOT->SetSelectedPad(0);
    TCanvas &c4 = d->AddTab("Quantiles");
    //c4.SetGrid();

    h2.SetFillColor(19);
    h1.SetFillColor(16);
    h0.SetFillColor(13);

    h2.SetTitle("Measured Residual-Distribution vs. Date");
    h2.SetYTitle("Residual [arcmin]");
    h2.SetStats(kFALSE);
    h2.GetXaxis()->SetTimeDisplay(1);
    h2.GetXaxis()->SetTimeFormat("%y/%m");
    h2.GetYaxis()->CenterTitle();
    h2.SetMinimum(0);
    h2.SetMaximum(h2.GetMaximum()*1.05);
    h2.SetMarkerStyle(kFullDotMedium);
    h2.GetXaxis()->SetRange(first, h3.GetNbinsX());
    h2.DrawCopy();

    //.DrawCopy("same");
    h1.DrawCopy("same");
    h0.DrawCopy("same");

    TGraph g;
    for (int i=0; i<h2.GetNbinsX(); i++)
    {
        Double_t x[2] = { h2.GetBinLowEdge(i+1), h2.GetBinLowEdge(i+1) };
        Double_t y[2] = { 0, h2.GetBinContent(i+1) };

        g.DrawGraph(2, x, y, "L");
    }

    DrawGrid(h2);

    gROOT->SetSelectedPad(0);
    TCanvas &c5 = d->AddTab("Combined");

    h2.DrawCopy();
    h1.DrawCopy("same");
    h0.DrawCopy("same");

    for (int i=0; i<h2.GetNbinsX(); i++)
    {
        Double_t x[2] = { h2.GetBinLowEdge(i+1), h2.GetBinLowEdge(i+1) };
        Double_t y[2] = { 0, h2.GetBinContent(i+1) };

        g.DrawGraph(2, x, y, "L");
    }

    g1.DrawClone("P");

    DrawGrid(h2);


}