#include #include #include #include #include #include #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]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[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 ×, Int_t i) { return TMath::BinarySearch(times.GetN(), times.GetX(), dates[i]);; } void DrawIndices(TGraph ×, 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("Q50", "", n-1, dates); TH1F h1("Q68", "", n-1, dates); TH1F h2("Q95", "", n-1, dates); TH1F h3("Q98", "", n-1, dates); Int_t first = 0; cout << "Number of bins: " << n << endl; cout << "Start Entries Median" << endl; MTime date; for (int i=0; dates[i+1]>0; i++) { Int_t bin0 = GetBin(g0, i); Int_t bin1 = GetBin(g0, i+1); if (bin0<0) bin0=0; if (bin0==bin1) continue; TMath::Quantiles(bin1-bin0, 4, g2.GetY()+bin0, quant, prob, kFALSE, NULL, 7); date.SetAxisTime(dates[i]); cout << date.GetStringFmt("%x", "de_DE") << " " << Form("%4d", bin1-bin0) << " " << Form("%4.1f", quant[0]*60) << endl; 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(); //h3.DrawCopy("same"); h1.DrawCopy("same"); h0.DrawCopy("same"); TGraph g; for (int i=0; iSetSelectedPad(0); TCanvas &c5 = d->AddTab("Combined"); h2.DrawCopy(); h1.DrawCopy("same"); h0.DrawCopy("same"); for (int i=0; i