Index: trunk/Mars/hawc/quality.C
===================================================================
--- trunk/Mars/hawc/quality.C	(revision 19921)
+++ trunk/Mars/hawc/quality.C	(revision 19921)
@@ -0,0 +1,405 @@
+#include <iostream>
+#include "TROOT.h"
+#include "TSystem.h"
+#include "TGraph.h"
+#include "TArrayD.h"
+#include "TCanvas.h"
+#include "TH1.h"
+#include "TAxis.h"
+#include "TText.h"
+#include "TLine.h"
+
+#include "fits.h"
+#include "MTime.h"
+
+#ifndef __CLING__
+#include <algorithm>
+#include <functional>
+#include <vector>
+#endif
+
+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(double start, TArrayD **vec, TString fname)
+{
+    fname += ".FTU_CONTROL_DAC.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    //cout << fname << endl;
+
+    Double_t time;
+    UShort_t th[4];
+
+    if (!file.SetPtrAddress("Time", &time))
+        return -1;
+
+    if (!file.SetPtrAddress("dac_ch", th))
+        return -1;
+
+    TGraph g;
+    g.SetName("Threshold");
+    while (file.GetNextRow())
+    {
+        if (Contains(vec, time, 10./(24*3600)))
+            g.SetPoint(g.GetN(), time*24*3600, th[0]);
+    }
+
+
+    gROOT->SetSelectedPad(0);
+    gPad->GetCanvas()->cd(1);
+    gPad->SetGrid();
+    gPad->SetTopMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    gPad->SetBottomMargin(0);
+
+    TH1C frame("Threshold_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
+    frame.SetMinimum(300);
+    frame.SetMaximum(1200);
+    frame.SetStats(kFALSE);
+    frame.GetXaxis()->SetTimeDisplay(true);
+    frame.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    frame.GetXaxis()->SetLabelSize(0.12);
+    frame.GetYaxis()->SetLabelSize(0.1);
+    frame.GetYaxis()->SetTitle("THRESHOLD");
+    frame.GetYaxis()->CenterTitle();
+    frame.GetYaxis()->SetTitleOffset(0.2);
+    frame.GetYaxis()->SetTitleSize(0.1);
+    frame.DrawCopy();
+
+    g.SetMarkerStyle(kFullDotMedium);
+    g.DrawClone("P");
+
+    return 0;
+}
+
+Int_t PlotCurrent(double start, TArrayD **vec, TString fname)
+{
+    fname += ".BIAS_CONTROL_AVERAGE_DATA.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    //cout << fname << endl;
+
+    Double_t time;
+    Float_t  Iavg[64];
+    Float_t  Tavg[64];
+    Float_t  Tpsu[2];
+    Float_t  Vavg[64];
+
+    if (!file.SetPtrAddress("Time", &time))
+        return -1;
+
+    if (!file.SetPtrAddress("Iavg", Iavg))
+        return -1;
+
+    if (!file.SetPtrAddress("Tavg", Tavg))
+        return -1;
+
+    if (!file.SetPtrAddress("Tavg_psu", Tpsu))
+        return -1;
+
+    if (!file.SetPtrAddress("Vavg", Vavg))
+        return -1;
+
+    TGraph g1;
+    TGraph g2;
+    TGraph g3;
+    TGraph g4;
+    g1.SetName("CurrentAvg");
+    g2.SetName("TemperatureAvg");
+    g3.SetName("TempPSUavg");
+    g4.SetName("VoltageAvg");
+
+    while (file.GetNextRow())
+        if (Contains(vec, time))
+        {
+            sort(Iavg, Iavg+64);
+            sort(Tavg, Tavg+64);
+            sort(Vavg, Vavg+64);
+
+            g1.SetPoint(g1.GetN(), time*24*3600, Iavg[31]);
+            g2.SetPoint(g2.GetN(), time*24*3600, Tavg[31]);
+            g3.SetPoint(g3.GetN(), time*24*3600, (Tpsu[0]+Tpsu[1])/2);
+            g4.SetPoint(g4.GetN(), time*24*3600, Vavg[31]);
+        }
+
+
+    gROOT->SetSelectedPad(0);
+    gPad->GetCanvas()->cd(2);
+    gPad->SetGrid();
+    gPad->SetTopMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    gPad->SetBottomMargin(0);
+
+    TH1C frame1("Current_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
+    frame1.SetMinimum(0);
+    frame1.SetMaximum(1000);
+    frame1.SetStats(kFALSE);
+    frame1.GetXaxis()->SetTimeDisplay(true);
+    frame1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    frame1.GetXaxis()->SetLabelSize(0.12);
+    frame1.GetYaxis()->SetLabelSize(0.1);
+    frame1.GetYaxis()->SetTitle("CURRENT");
+    frame1.GetYaxis()->CenterTitle();
+    frame1.GetYaxis()->SetTitleOffset(0.2);
+    frame1.GetYaxis()->SetTitleSize(0.1);
+    frame1.DrawCopy();
+
+    g1.SetMarkerStyle(kFullDotMedium);
+    g1.DrawClone("P");
+
+
+
+    gROOT->SetSelectedPad(0);
+    gPad->GetCanvas()->cd(5);
+    gPad->SetGrid();
+    gPad->SetTopMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    gPad->SetBottomMargin(0);
+
+    TH1C frame4("Voltage_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
+    frame4.SetMinimum(29);
+    frame4.SetMaximum(30);
+    frame4.SetStats(kFALSE);
+    frame4.GetXaxis()->SetTimeDisplay(true);
+    frame4.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    frame4.GetXaxis()->SetLabelSize(0.12);
+    frame4.GetYaxis()->SetLabelSize(0.1);
+    frame4.GetYaxis()->SetTitle("VOLTAGE");
+    frame4.GetYaxis()->CenterTitle();
+    frame4.GetYaxis()->SetTitleOffset(0.2);
+    frame4.GetYaxis()->SetTitleSize(0.1);
+    frame4.DrawCopy();
+
+    g4.SetMarkerStyle(kFullDotMedium);
+    g4.DrawClone("P");
+
+
+
+    gROOT->SetSelectedPad(0);
+    gPad->GetCanvas()->cd(6);
+    gPad->SetGrid();
+    gPad->SetTopMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+    gPad->SetBottomMargin(0);
+
+    TH1C frame2("Temperature_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
+    frame2.SetMinimum(0);
+    frame2.SetMaximum(40);
+    frame2.SetStats(kFALSE);
+    frame2.GetXaxis()->SetTimeDisplay(true);
+    frame2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    frame2.GetXaxis()->SetLabelSize(0.12);
+    frame2.GetYaxis()->SetLabelSize(0.1);
+    frame2.GetYaxis()->SetTitle("TEMPERATURE");
+    frame2.GetYaxis()->CenterTitle();
+    frame2.GetYaxis()->SetTitleOffset(0.2);
+    frame2.GetYaxis()->SetTitleSize(0.1);
+    frame2.DrawCopy();
+
+    g2.SetMarkerStyle(kFullDotMedium);
+    g2.DrawClone("P");
+
+    g3.SetMarkerColor(kBlue);
+    g3.SetMarkerStyle(kFullDotMedium);
+    g3.DrawClone("P");
+
+    return 0;
+}
+
+Int_t PlotRate(double start, TArrayD **vec, TString fname)
+{
+    fname += ".FTM_CONTROL_DATA.fits";
+
+    fits file(fname.Data());
+    if (!file)
+    {
+        cerr << fname << ": " << gSystem->GetError() << endl;
+        return -2;
+    }
+
+    //cout << fname << endl;
+
+    Double_t time;
+    uint32_t trg_counter;
+    uint32_t dead_time, run_time;
+
+    if (!file.SetPtrAddress("Time",  &time))
+        return -1;
+
+    if (!file.SetPtrAddress("trg_counter", &trg_counter))
+        return -1;
+    if (!file.SetPtrAddress("dead_time", &dead_time))
+        return -1;
+    if (!file.SetPtrAddress("run_time", &run_time))
+        return -1;
+
+    TGraph g1, g2;
+    g1.SetName("TriggerRate");
+    g2.SetName("RelOnTime");
+
+    uint32_t prev = 0;
+
+    while (file.GetNextRow())
+        if (Contains(vec, time))
+        {
+            g1.SetPoint(g1.GetN(), time*24*3600, 1e8*(trg_counter-prev)/(run_time-dead_time));
+            g2.SetPoint(g2.GetN(), time*24*3600, 1.-float(dead_time)/run_time);
+            prev = trg_counter;
+        }
+
+
+    gROOT->SetSelectedPad(0);
+    gPad->GetCanvas()->cd(3);
+    gPad->SetGrid();
+    gPad->SetTopMargin(0);
+    gPad->SetBottomMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+
+    TH1C frame1("Rate_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
+    frame1.SetMinimum(0);
+    frame1.SetMaximum(90);
+    frame1.SetStats(kFALSE);
+    frame1.GetXaxis()->SetTimeDisplay(true);
+    frame1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    frame1.GetXaxis()->SetLabelSize(0.12);
+    frame1.GetYaxis()->SetLabelSize(0.1);
+    frame1.GetYaxis()->SetTitle("TRIGGER RATE");
+    frame1.GetYaxis()->CenterTitle();
+    frame1.GetYaxis()->SetTitleOffset(0.2);
+    frame1.GetYaxis()->SetTitleSize(0.1);
+    frame1.DrawCopy();
+
+    g1.SetMarkerStyle(kFullDotMedium);
+    g1.DrawClone("P");
+
+
+
+    gROOT->SetSelectedPad(0);
+    gPad->GetCanvas()->cd(4);
+    gPad->SetGrid();
+    gPad->SetTopMargin(0);
+    gPad->SetBottomMargin(0);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+
+    TH1C frame2("OnTime_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
+    frame2.SetMinimum(0);
+    frame2.SetMaximum(1);
+    frame2.SetStats(kFALSE);
+    frame2.GetXaxis()->SetTimeDisplay(true);
+    frame2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    frame2.GetXaxis()->SetLabelSize(0.12);
+    frame2.GetYaxis()->SetLabelSize(0.1);
+    frame2.GetYaxis()->SetTitle("RELATIVE ON TIME");
+    frame2.GetYaxis()->CenterTitle();
+    frame2.GetYaxis()->SetTitleOffset(0.2);
+    frame2.GetYaxis()->SetTitleSize(0.1);
+    frame2.DrawCopy();
+
+    g2.SetMarkerStyle(kFullDotMedium);
+    g2.DrawClone("P");
+
+    return 0;
+}
+
+int quality(const char *basepath, UInt_t y=0, UInt_t m=0, UInt_t d=0, const char *outpath=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("%s/%04d/%02d/%02d/%04d%02d%02d", basepath, y, m, d, y, m, d);
+
+    cout << "quality" << endl;
+    cout << "-------" << endl;
+    cout << endl;
+    cout << "Night: " << Form("%04d-%02d-%02d", y, m, d) << endl;
+    cout << endl;
+
+    double start = MTime(y, m, d).GetMjd()-40587;
+
+    TArrayD run, beg, end;
+
+    TArrayD *runs[3] = { &run, &beg, &end };
+
+    //check if the sqm was already installed on the telescope                                                                                                       
+    TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 1120);
+    c->Divide(1, 7, 1e-5, 1e-5);
+
+    PlotThresholds(start, runs, fname);
+    PlotCurrent(start, runs, fname);
+    PlotRate(start, runs, fname);
+
+    gROOT->SetSelectedPad(0);
+    gPad->GetCanvas()->cd(7);
+    gPad->SetTopMargin(0);
+    gPad->SetBottomMargin(0.99);
+    gPad->SetRightMargin(0.001);
+    gPad->SetLeftMargin(0.04);
+
+    TH1C frame1("Axis", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
+    frame1.SetMinimum(0);
+    frame1.SetMaximum(1);
+    frame1.SetStats(kFALSE);
+    frame1.GetXaxis()->SetTimeDisplay(true);
+    frame1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 17:00:00 GMT");
+    frame1.GetXaxis()->SetLabelSize(0.2);
+    frame1.GetXaxis()->SetTitleSize(0.2);
+    frame1.GetXaxis()->SetTitle("HAWC Local Time");
+    frame1.GetXaxis()->CenterTitle();
+    frame1.DrawCopy();
+
+    if (outpath)
+        c->SaveAs(outpath);
+
+    return 0;
+}