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

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