Index: trunk/Mars/fact/plots/quality.C
===================================================================
--- trunk/Mars/fact/plots/quality.C	(revision 15244)
+++ trunk/Mars/fact/plots/quality.C	(revision 15244)
@@ -0,0 +1,491 @@
+Bool_t Contains(TArrayD **vec, Double_t t0, Double_t range=0)
+{
+    TArrayD *arr0 = vec[0];
+    TArrayD *arr1 = vec[1];
+    TArrayD *arr2 = vec[2];
+
+    for (int i=0; i<arr0->GetSize(); i++)
+    {
+        Double_t start = (*arr1)[i];
+        Double_t stop  = (*arr2)[i];
+
+        if (stop>start+305./24/3600)
+            stop = start+305./24/3600;
+
+        if (t0>start-range && t0<stop+range)
+        //{
+        //    if (fmod(t0, 1)>4./24 && fmod(t0,1)<4.1/24)
+        //        cout << t0-start << " " <<t0 << " " << stop-t0 << " " << start-15779 << " " << stop-15779 << " " << (*arr0)[i] << endl;
+            return kTRUE;
+        //}
+    }
+
+    return arr0->GetSize()==0;
+}
+
+Int_t PlotThresholds(TArrayD **vec, TString fname)
+{
+    fname += ".RATE_CONTROL_THRESHOLD.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    //cout << fname << endl;
+
+    Double_t time;
+    UShort_t th;
+
+    if (!file.SetPtrAddress("Time", &time))
+        return -1;
+
+    if (!file.SetPtrAddress("threshold", &th))
+        return -1;
+
+    TGraph g;
+
+    while (file.GetNextRow())
+    {
+        if (Contains(vec, time, 10./(24*3600)))
+            g.SetPoint(g.GetN(), time*24*3600, th);
+    }
+
+    g.SetMarkerStyle(kFullDotMedium);
+    g.GetXaxis()->SetTimeDisplay(true);
+    g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    g.GetXaxis()->SetLabelSize(0.12);
+    g.GetYaxis()->SetLabelSize(0.1);
+    g.GetYaxis()->SetTitle("THRESHOLD");
+    g.GetYaxis()->SetTitleOffset(0.2);
+    g.GetYaxis()->SetTitleSize(0.1);
+    g.DrawClone("AP");
+
+    return 0;
+}
+
+#include <vector>
+#include <pair>
+
+vector<pair<double, ln_equ_posn>> vec;
+
+ln_equ_posn FindPointing(Double_t time)
+{
+    for (int i=0; i<vec.size(); i++)
+        if (time<vec[i].first)
+            return i==0 ? ln_equ_posn() : vec[i-1].second;
+
+    return vec[vec.size()-1].second;
+}
+
+Float_t Prediction(Double_t time)
+{
+    Double_t jd = time + 40587 + 2400000.5;
+
+    Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
+
+    // Nova::EquPosn pos = FindPointing(time);
+
+    // get local position of moon
+    Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
+
+    // Distance between source and moon
+    //const double angle =
+    //    MAstro::AngularDistance(TMath::DegToRad()*(90-moon.dec),
+    //                            TMath::DegToRad()*moon.ra,
+    //                            TMath::DegToRad()*(90-pos.dec),
+    //                            TMath::DegToRad()*pos.ra);
+
+    // Current prediction
+    //double cang = sin(angle);
+    double calt = sin(hrzm.alt*TMath::DegToRad());
+
+    double disk = Nova::GetLunarDisk(jd);
+
+    double lc = /*cang==0 ? 0 :*/ calt*pow(disk, 2.5);///sqrt(cang);
+
+    return 5+103*lc>4.5 ? 5+103*lc : 4.5;
+}
+
+Int_t ReadSources(TString fname)
+{
+    fname += ".DRIVE_CONTROL_SOURCE_POSITION.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    Double_t time, ra, dec;
+    if (!file.SetPtrAddress("Time",  &time))
+        return -1;
+    if (!file.SetPtrAddress("Ra_cmd", &ra))
+        return -1;
+    if (!file.SetPtrAddress("Dec_cmd", &dec))
+        return -1;
+
+    while (file.GetNextRow())
+    {
+        ln_equ_posn p;
+        p.ra  = ra*15;
+        p.dec = dec;
+
+        vec.push_back(make_pair(time, p));
+    }
+
+    return 0;
+}
+
+Int_t PlotCurrent(TArrayD **vec, TString fname)
+{
+    Int_t rc = ReadSources(fname);
+    if (rc<0)
+        return rc;
+
+    fname += ".FEEDBACK_CALIBRATED_CURRENTS.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    //cout << fname << endl;
+
+    Double_t time;
+    Float_t  Imed;
+
+    if (!file.SetPtrAddress("Time",  &time))
+        return -1;
+
+    if (!file.SetPtrAddress("I_med", &Imed))
+        return -1;
+
+    TGraph g1;
+    TGraph g2;
+
+    while (file.GetNextRow())
+        if (Contains(vec, time))
+        {
+            g1.SetPoint(g1.GetN(), time*24*3600, Imed);
+            g2.SetPoint(g2.GetN(), time*24*3600, Prediction(time));
+        }
+
+    g1.SetMinimum(0);
+    g1.SetMaximum(99);
+
+    g1.SetMarkerStyle(kFullDotMedium);
+    g1.GetXaxis()->SetTimeDisplay(true);
+    g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    g1.GetXaxis()->SetLabelSize(0.12);
+    g1.GetYaxis()->SetLabelSize(0.1);
+    g1.GetYaxis()->SetTitle("CURRENT");
+    g1.GetYaxis()->SetTitleOffset(0.2);
+    g1.GetYaxis()->SetTitleSize(0.1);
+    g1.DrawClone("AP");
+
+    g2.SetMarkerColor(kBlue);
+    g2.SetMarkerStyle(kFullDotMedium);
+    g2.GetXaxis()->SetTimeDisplay(true);
+    g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    g2.GetXaxis()->SetLabelSize(0.12);
+    g2.GetYaxis()->SetLabelSize(0.1);
+    g2.GetYaxis()->SetTitle("CURRENT");
+    g2.GetYaxis()->SetTitleOffset(0.2);
+    g2.GetYaxis()->SetTitleSize(0.1);
+    g2.DrawClone("P");
+
+    g1.DrawClone("P");
+
+    return 0;
+}
+
+Int_t PlotRate(TArrayD **vec, TString fname)
+{
+    fname += ".FTM_CONTROL_TRIGGER_RATES.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    //cout << fname << endl;
+
+    Double_t time;
+    Float_t  rate;
+    Float_t  ontime, elapsed;
+
+    if (!file.SetPtrAddress("Time",  &time))
+        return -1;
+
+    if (!file.SetPtrAddress("TriggerRate", &rate))
+        return -1;
+    if (!file.SetPtrAddress("OnTime", &ontime))
+        return -1;
+    if (!file.SetPtrAddress("ElapsedTime", &elapsed))
+        return -1;
+
+    TGraph g1, g2;
+
+    while (file.GetNextRow())
+        if (Contains(vec, time))
+        {
+            if (rate>=0)
+            {
+                g1.SetPoint(g1.GetN(), time*24*3600, rate);
+                g2.SetPoint(g2.GetN(), time*24*3600, ontime/elapsed);
+            }
+        }
+
+    g1.SetMinimum(0);
+    g1.SetMaximum(200);
+    g1.SetMarkerStyle(kFullDotMedium);
+    g1.GetXaxis()->SetTimeDisplay(true);
+    g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    g1.GetXaxis()->SetLabelSize(0.12);
+    g1.GetYaxis()->SetLabelSize(0.1);
+    g1.GetYaxis()->SetTitle("TRIGGER RATE");
+    g1.GetYaxis()->SetTitleOffset(0.2);
+    g1.GetYaxis()->SetTitleSize(0.1);
+    g1.DrawClone("AP");
+
+    gROOT->SetSelectedPad(0);
+    gPad->GetCanvas()->cd(4);
+
+    gPad->SetGridx();
+    gPad->SetTopMargin(0);
+    gPad->SetBottomMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+
+    g2.SetMinimum(0);
+    g2.SetMaximum(1);
+    g2.SetMarkerStyle(kFullDotMedium);
+    g2.GetXaxis()->SetTimeDisplay(true);
+    g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    g2.GetXaxis()->SetLabelSize(0.12);
+    g2.GetYaxis()->SetLabelSize(0.1);
+    g2.GetYaxis()->SetTitle("RELATIVE ON TIME");
+    g2.GetYaxis()->SetTitleOffset(0.2);
+    g2.GetYaxis()->SetTitleSize(0.1);
+    g2.DrawClone("AP");
+
+
+    return 0;
+}
+
+Int_t PlotPointing(TArrayD **vec, TString fname)
+{
+    fname += ".DRIVE_CONTROL_POINTING_POSITION.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    //cout << fname << endl;
+
+    Double_t time;
+    Double_t zd;
+
+    if (!file.SetPtrAddress("Time",  &time))
+        return -1;
+    if (!file.SetPtrAddress("Zd", &zd))
+        return -1;
+
+    TGraph g;
+
+    while (file.GetNextRow())
+        if (Contains(vec, time))
+            g.SetPoint(g.GetN(), time*24*3600, 90-zd);
+
+    g.SetMinimum(0);
+    g.SetMaximum(90);
+    g.SetMarkerStyle(kFullDotMedium);
+    g.GetXaxis()->SetTimeDisplay(true);
+    g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    g.GetXaxis()->SetLabelSize(0.1);
+    g.GetYaxis()->SetLabelSize(0.1);
+    g.GetYaxis()->SetTitle("ELEVATION");
+    g.GetYaxis()->SetTitleOffset(0.2);
+    g.GetYaxis()->SetTitleSize(0.1);
+    g.DrawClone("AP");
+
+    return 0;
+}
+/*
+Int_t PlotMem(TArrayD **vec, TString fname)
+{
+    fname += ".FAD_CONTROL_STATISTICS1.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    //cout << fname << endl;
+
+    Double_t time;
+    UInt_t buf[5];
+    ULong64_t mem[4];
+    UInt_t rate[2];
+
+    if (!file.SetPtrAddress("Time",  &time))
+        return -1;
+    if (!file.SetPtrAddress("bufferInfo", buf) ||
+        !file.SetPtrAddress("memInfo",    mem) ||
+        !file.SetPtrAddress("rateNew",    rate))
+        return -1;
+
+    TGraph g;
+
+    while (file.GetNextRow())
+        if (Contains(vec, time))
+            g.SetPoint(g.GetN(), time*24*3600, mem[1]);
+
+    g.SetMarkerStyle(kFullDotMedium);
+    g.GetXaxis()->SetTimeDisplay(true);
+    g.GetXaxis()->SetTimeFormat("%H:%M");
+    g.GetXaxis()->SetTimeOffset(0, "gmt");
+    g.GetXaxis()->SetLabelSize(0.12);
+    g.GetYaxis()->SetLabelSize(0.1);
+    //g.GetYaxis()->SetTitle("ELEVATION");
+    g.GetYaxis()->SetTitleOffset(0.2);
+    g.GetYaxis()->SetTitleSize(0.1);
+    g.DrawClone("AP");
+
+    return 0;
+}
+*/
+void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0)
+{
+    // To get correct dates in the histogram you have to add
+    // the MJDREF offset (should be 40587) and 9131.
+
+    if (y==0)
+    {
+        UInt_t nt = MTime(MTime(-1).GetMjd()-1.5).GetNightAsInt();
+        y = nt/10000;
+        m = (nt/100)%100;
+        d = nt%100;
+
+        cout << y << "/" << m << "/" << d << endl;
+    }
+
+    TString fname=Form("/fact/aux/%04d/%02d/%02d/%04d%02d%02d", y, m, d, y, m, d);
+
+    UInt_t night = MTime(y, m, d, 0).GetNightAsInt();
+
+    MSQLMagic serv("sql.rc");
+    Bool_t con = serv.IsConnected();
+
+    cout << "quality" << endl;
+    cout << "-------" << endl;
+    cout << endl;
+    if (con)
+    {
+        cout << "Connected to " << serv.GetName() << endl;
+        cout << endl;
+    }
+    cout << "Night: " << Form("%04d-%02d-%02d", y, m, d) << endl;
+    cout << endl;
+
+    TArrayD run, beg, end;
+
+    TArrayD *runs[3] = { &run, &beg, &end };
+
+    if (con)
+    {
+        TString query;
+        query += "SELECT fRunID, fRunStart, fRunStop FROM runinfo";
+        query += " WHERE fNight=";
+        query += night;
+        query += " AND fRunTypeKEY=1 ORDER BY fRunStart";
+
+        TSQLResult *res = serv.Query(query);
+        if (!res)
+            return kFALSE;
+
+        run.Set(res->GetRowCount());
+        beg.Set(res->GetRowCount());
+        end.Set(res->GetRowCount());
+
+        Int_t n = 0;
+
+        TSQLRow *row = 0;
+        while ((row=res->Next()))
+        {
+            run[n] = atoi((*row)[0]);
+            beg[n] = MTime((*row)[1]).GetMjd()-40587;
+            end[n] = MTime((*row)[2]).GetMjd()-40587;
+            n++;
+            delete row;
+        }
+
+        delete res;
+
+        if (n==0)
+            cout << "WARNING - No data runs in db, displaying all data." << endl;
+    }
+
+    TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 960);
+    c->Divide(1, 5, 1e-5, 1e-5);
+
+    gROOT->SetSelectedPad(0);
+    c->cd(1);
+    gPad->SetGridx();
+    gPad->SetTopMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    gPad->SetBottomMargin(0);
+    cout << PlotThresholds(runs, fname) << endl;
+
+    gROOT->SetSelectedPad(0);
+    c->cd(2);
+    gPad->SetGridx();
+    gPad->SetTopMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    gPad->SetBottomMargin(0);
+    cout << PlotCurrent(runs, fname) << endl;
+
+    gROOT->SetSelectedPad(0);
+    c->cd(3);
+    gPad->SetGridx();
+    gPad->SetTopMargin(0);
+    gPad->SetBottomMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    cout << PlotRate(runs, fname) << endl;
+
+    gROOT->SetSelectedPad(0);
+    c->cd(5);
+    gPad->SetGridx();
+    gPad->SetTopMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    cout << PlotPointing(runs, fname) << endl;
+/*
+    gROOT->SetSelectedPad(0);
+    c->cd(6);
+    gPad->SetGridx();
+    gPad->SetTopMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    cout << PlotMem(runs, fname) << endl;
+*/
+    c->SaveAs(Form("%04d%02d%02d-quality.pdf", y, m, d));
+}
+
+//20130314_141