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

Last change on this file since 10088 was 10088, checked in by tbretz, 10 years ago
Added output of filename for outliers.
File size: 6.5 KB
Line 
1#include <iostream>
2#include <iomanip>
3#include <fstream>
4
5#include <TGraph.h>
6#include <TCanvas.h>
7#include <TSystem.h>
8
9#include "mars/MTime.h"
10#include "mars/MMath.h"
11#include "mars/MAstro.h"
12#include "mars/MDirIter.h"
13#include "mars/MStatusDisplay.h"
14
15using namespace std;
16
17int ReadFile(const char *fname, TGraph &g1, TGraph &g2, TGraph &g0)
18{
19    ifstream fin(fname);
20
21    cout << "Reading " << setw(23) << gSystem->BaseName(fname) << "..." << flush;
22
23    MTime t;
24    while (1)
25    {
26        TString str;
27        str.ReadLine(fin);
28        if (!fin)
29            break;
30
31        if (str.Contains("#"))
32            continue;
33
34        Double_t alt, az, dalt, daz, mjd;
35        if (sscanf(str.Data(), "%lf %lf %*f %*f %*f %*f %lf %lf %lf",
36                   &az, &alt, &dalt, &daz, &mjd)!=5)
37            continue;
38
39        if (dalt==0)
40            continue;
41
42        Double_t res = MAstro::GetDevAbs(90-(alt-dalt), -dalt, -daz)*60;
43
44        if (res>12)
45            continue;
46
47        t.SetMjd(mjd);
48
49        g0.SetPoint(g0.GetN(), t.GetAxisTime(), g0.GetN());
50        g1.SetPoint(g1.GetN(), t.GetAxisTime(), res);
51        g2.SetPoint(g2.GetN(), g2.GetN(),       res);
52
53        if (g0.GetN()<2 || g0.GetX()[g0.GetN()-2]<g0.GetX()[g0.GetN()-1])
54            continue;
55
56        cout << "WARNING-TIME-SKEW! " << mjd << " " << flush;
57    }
58
59    cout << g1.GetN() << endl;
60
61    return g1.GetN();
62}
63
64void DrawGrid(TH1 &h)
65{
66    TGraph l;
67    l.SetLineStyle(kDashed);
68    l.SetLineWidth(2);
69    l.SetLineColor(kBlue);
70
71    Int_t style[] = { kSolid, 9, 7 };
72
73    Double_t x[2] = { h.GetXaxis()->GetXmin(), h.GetXaxis()->GetXmax() };
74    Double_t y[2];
75
76    for (int i=0; i<3; i++)
77    {
78        y[0] = y[1] = i+1;
79        l.SetLineStyle(style[i]);
80        l.DrawGraph(2, x, y, "L");
81    }
82
83    if (x[0]<1e6)
84        return;
85
86    Double_t X[2];
87    Double_t Y[2] = { h.GetMinimum(), h.GetMaximum() };
88
89    l.SetLineWidth(1);
90
91    for (int i=2004; i<2020; i++)
92    {
93        for (int j=0; j<12; j++)
94        {
95            X[0] = X[1] = MTime(i, 1+j, 1).GetAxisTime();
96
97            if (X[0]<x[0] || X[0]>x[1])
98                continue;
99
100            l.SetLineStyle(j==0 ? kDashed : kDotted);
101            l.DrawGraph(2, X, Y, "L");
102        }
103    }
104}
105
106void DrawTimes(TH1 &h)
107{
108    TGraph l;
109    l.SetLineColor(kBlue);
110
111    Double_t x[2];
112    Double_t y[2] = { h.GetMinimum(), h.GetMaximum() };
113
114    for (int i=0; dates[i+1]>0; i++)
115    {
116        x[0] = x[1] = dates[i];
117        l.DrawGraph(2, x, y, "L");
118    }
119
120    DrawGrid(h);
121}
122
123Int_t GetBin(TGraph &times, Int_t i)
124{
125    return TMath::BinarySearch(times.GetN(), times.GetX(), dates[i]);;
126}
127
128void DrawIndices(TGraph &times, TH1 &h)
129{
130    TGraph l;
131    l.SetLineStyle(kDashed);
132    l.SetLineColor(kBlue);
133
134    Double_t x[2];
135    Double_t y[2] = { h.GetMinimum(), h.GetMaximum() };
136
137    for (int i=0; dates[i+1]>0; i++)
138    {
139        Int_t bin = GetBin(times, i);
140        if (bin<0)
141            continue;
142
143        x[0] = x[1] = bin+0.5;
144        l.DrawGraph(2, x, y, "L");
145    }
146
147    DrawGrid(h);
148}
149
150void plot(MDirIter &Next)
151{
152    Next.Sort();
153
154    TGraph g1;
155    TGraph g2;
156    TGraph g0;
157
158    while (1)
159    {
160        TString name=Next();
161        if (name.IsNull())
162            break;
163
164        ReadFile(name, g1, g2, g0);
165    }
166
167    MStatusDisplay *d = new MStatusDisplay;
168
169    gROOT->SetSelectedPad(0);
170    TCanvas &c0 = d->AddTab("TimeLine");
171    c0.SetGrid();
172
173    g0.SetTitle("TimeLine");
174    g0.GetXaxis()->SetTimeDisplay(1);
175    g0.SetMarkerStyle(kFullDotMedium);
176    g0.SetLineWidth(1);
177    g0.SetLineColor(kGreen);
178    g0.DrawClone("ALP");
179
180    DrawTimes(*g0.GetHistogram());
181
182
183    gROOT->SetSelectedPad(0);
184    TCanvas &c1 = d->AddTab("ResDate");
185    //c1.SetGrid();
186
187    g1.SetTitle("Measured Residual vs. Date");
188    g1.GetXaxis()->SetTimeDisplay(1);
189//    g1.SetMaximum(10);
190    g1.SetMarkerStyle(kFullDotMedium);
191    g1.DrawClone("AP");
192
193    DrawTimes(*g1.GetHistogram());
194
195
196    gROOT->SetSelectedPad(0);
197    TCanvas &c2 = d->AddTab("ResIdx");
198    c2.SetGrid();
199
200    g2.SetTitle("Measured Residual vs. Index");
201//    g2.SetMaximum(0.1);
202    g2.SetMarkerStyle(kFullDotMedium);
203    g2.DrawClone("AP");
204
205    DrawIndices(g0, *g2.GetHistogram());
206
207    Double_t prob[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
208    Double_t quant[4];
209
210    int n;
211    for (n=0; dates[n]>0; n++);
212
213    TH1F h0("Q50", "", n-1, dates);
214    TH1F h1("Q68", "", n-1, dates);
215    TH1F h2("Q95", "", n-1, dates);
216    TH1F h3("Q98", "", n-1, dates);
217
218    Int_t first = 0;
219
220    cout << "Number of bins: " << n << endl;
221    cout << "Start      Entries Median" << endl;
222
223    MTime date;
224    for (int i=0; dates[i+1]>0; i++)
225    {
226        Int_t bin0 = GetBin(g0, i);
227        Int_t bin1 = GetBin(g0, i+1);
228
229        //cout << i << " " << bin0 << " " << bin1 << endl;
230        if (bin0<0)
231            continue;
232
233        if (bin0==bin1)
234            continue;
235
236        TMath::Quantiles(bin1-bin0, 4, g2.GetY()+bin0, quant, prob, kFALSE, NULL, 7);
237
238        date.SetAxisTime(dates[i]);
239        cout << date.GetStringFmt("%x", "de_DE") << "  " << Form("%4d", bin1-bin0) << "   " << Form("%4.1f", quant[0]*60) << endl;
240
241        h0.SetBinContent(i+1, quant[0]);
242        h1.SetBinContent(i+1, quant[1]);
243        h2.SetBinContent(i+1, quant[2]);
244        h3.SetBinContent(i+1, quant[3]);
245
246        if (first==0)
247            first = i;
248    }
249
250
251    gROOT->SetSelectedPad(0);
252    TCanvas &c4 = d->AddTab("Quantiles");
253    //c4.SetGrid();
254
255    h2.SetFillColor(19);
256    h1.SetFillColor(16);
257    h0.SetFillColor(13);
258
259    h2.SetTitle("Measured Residual-Distribution vs. Date");
260    h2.SetYTitle("Residual [arcmin]");
261    h2.SetStats(kFALSE);
262    h2.GetXaxis()->SetTimeDisplay(1);
263    h2.GetXaxis()->SetTimeFormat("%y/%m");
264    h2.GetYaxis()->CenterTitle();
265    h2.SetMinimum(0);
266    h2.SetMaximum(h2.GetMaximum()*1.05);
267    h2.SetMarkerStyle(kFullDotMedium);
268    h2.GetXaxis()->SetRange(first, h3.GetNbinsX());
269    h2.DrawCopy();
270
271    //h3.DrawCopy("same");
272    h1.DrawCopy("same");
273    h0.DrawCopy("same");
274
275    TGraph g;
276    for (int i=0; i<h2.GetNbinsX(); i++)
277    {
278        Double_t x[2] = { h2.GetBinLowEdge(i+1), h2.GetBinLowEdge(i+1) };
279        Double_t y[2] = { 0, h2.GetBinContent(i+1) };
280
281        g.DrawGraph(2, x, y, "L");
282    }
283
284    DrawGrid(h2);
285
286    gROOT->SetSelectedPad(0);
287    TCanvas &c5 = d->AddTab("Combined");
288
289    h2.DrawCopy();
290    h1.DrawCopy("same");
291    h0.DrawCopy("same");
292
293    for (int i=0; i<h2.GetNbinsX(); i++)
294    {
295        Double_t x[2] = { h2.GetBinLowEdge(i+1), h2.GetBinLowEdge(i+1) };
296        Double_t y[2] = { 0, h2.GetBinContent(i+1) };
297
298        g.DrawGraph(2, x, y, "L");
299    }
300
301    g1.DrawClone("P");
302
303    DrawGrid(h2);
304}
Note: See TracBrowser for help on using the repository browser.