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 


15  using namespace std;


16 


17  int 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(altdalt), 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 << "WARNINGTIMESKEW! " << mjd << " " << flush;


57  }


58 


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


60 


61  return g1.GetN();


62  }


63 


64  void 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 


106  void 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 


123  Int_t GetBin(TGraph ×, Int_t i)


124  {


125  return TMath::BinarySearch(times.GetN(), times.GetX(), dates[i]);;


126  }


127 


128  void DrawIndices(TGraph ×, 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 


150  void 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("Q98", "", n1, dates);


214  TH1F h1("Q95", "", n1, dates);


215  TH1F h2("Q68", "", n1, dates);


216  TH1F h3("Q50", "", n1, dates);


217 


218  Int_t first = 0;


219 


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


221 


222  for (int i=0; dates[i+1]>0; i++)


223  {


224  Int_t bin0 = GetBin(g0, i);


225  Int_t bin1 = GetBin(g0, i+1);


226 


227  cout << i << " " << bin0 << " " << bin1 << endl;


228  if (bin0<0)


229  continue;


230 


231  if (bin0==bin1)


232  continue;


233 


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


235 


236  h0.SetBinContent(i+1, quant[0]);


237  h1.SetBinContent(i+1, quant[1]);


238  h2.SetBinContent(i+1, quant[2]);


239  h3.SetBinContent(i+1, quant[3]);


240 


241  if (first==0)


242  first = i;


243  }


244 


245 


246  gROOT>SetSelectedPad(0);


247  TCanvas &c4 = d>AddTab("Quantiles");


248  //c4.SetGrid();


249 


250  h2.SetFillColor(19);


251  h1.SetFillColor(16);


252  h0.SetFillColor(13);


253 


254  h2.SetTitle("Measured ResidualDistribution vs. Date");


255  h2.SetYTitle("Residual [arcmin]");


256  h2.SetStats(kFALSE);


257  h2.GetXaxis()>SetTimeDisplay(1);


258  h2.GetXaxis()>SetTimeFormat("%y/%m");


259  h2.GetYaxis()>CenterTitle();


260  h2.SetMinimum(0);


261  h2.SetMaximum(h2.GetMaximum()*1.05);


262  h2.SetMarkerStyle(kFullDotMedium);


263  h2.GetXaxis()>SetRange(first, h3.GetNbinsX());


264  h2.DrawCopy();


265 


266  //.DrawCopy("same");


267  h1.DrawCopy("same");


268  h0.DrawCopy("same");


269 


270  TGraph g;


271  for (int i=0; i<h2.GetNbinsX(); i++)


272  {


273  Double_t x[2] = { h2.GetBinLowEdge(i+1), h2.GetBinLowEdge(i+1) };


274  Double_t y[2] = { 0, h2.GetBinContent(i+1) };


275 


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


277  }


278 


279  DrawGrid(h2);


280 


281  gROOT>SetSelectedPad(0);


282  TCanvas &c5 = d>AddTab("Combined");


283 


284  h2.DrawCopy();


285  h1.DrawCopy("same");


286  h0.DrawCopy("same");


287 


288  for (int i=0; i<h2.GetNbinsX(); i++)


289  {


290  Double_t x[2] = { h2.GetBinLowEdge(i+1), h2.GetBinLowEdge(i+1) };


291  Double_t y[2] = { 0, h2.GetBinContent(i+1) };


292 


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


294  }


295 


296  g1.DrawClone("P");


297 


298  DrawGrid(h2);


299 


300 


301  }

