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

Last change on this file since 12582 was 10088, checked in by tbretz, 14 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.