Index: trunk/Mars/fact/analysis/gain/extract_aux.C
===================================================================
--- trunk/Mars/fact/analysis/gain/extract_aux.C	(revision 17370)
+++ trunk/Mars/fact/analysis/gain/extract_aux.C	(revision 17370)
@@ -0,0 +1,377 @@
+#include <iostream>
+
+#include <TString.h>
+#include <TFile.h>
+#include <TGraph.h>
+
+#include "zfits.h"
+#include "PixelMap.h"
+#include "Interpolator2D.h"
+
+#include "MMath.h"
+#include "MHCamera.h"
+#include "MGeomCamFACT.h"
+
+PixelMap pmap;
+
+bool extract_current(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
+{
+    TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FEEDBACK_CALIBRATED_CURRENTS.fits",
+                        day/10000, (day/100)%100, day%100, day);
+
+    // =====================================================================
+
+    // Now we read the temperatures of the sensors during the
+    // single pe run and average them
+    fits faux(naux.Data());
+    if (!faux)
+    {
+        cout << "Could not open " << naux << endl;
+        return false;
+    }
+
+    Double_t time;
+    Float_t  I[416];
+
+    if (!faux.SetPtrAddress("Time",  &time))
+        return false;
+    if (!faux.SetPtrAddress("I", I))
+        return false;
+
+    TArrayD avg(320);
+    TArrayD rms(320);
+    TArrayI cnt(320);
+
+    while (faux.GetNextRow())
+    {
+        if (time<tstart || time>tstop)
+            continue;
+
+        TArrayD med(320);
+
+        double mean = 0;
+
+        int j=0;
+        for (int i=0; i<320; i++)
+        {
+            if (i!=66 && i!=191 && i!=193 && i!=17 && i!=206)
+                med[j++]  = I[i];
+
+            avg[i] += I[i];
+            rms[i] += I[i]*I[i];
+            cnt[i]++;
+        }
+
+        Double_t median;
+        MMath::MedianDev(315, med.GetArray(), median);
+        g.SetPoint(g.GetN(), time, median);
+    }
+
+    if (cnt[0]==0)
+        return 2;
+
+    for (int i=0; i<320; i++)
+    {
+        avg[i] /= cnt[i];
+        rms[i] /= cnt[i];
+        rms[i] -= avg[i]*avg[i];
+        rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
+    }
+
+    double mn = 0;
+    for (int i=0; i<1440; i++)
+    {
+        mn += avg[pmap.index(i).hv()];
+        havg.SetBinContent(i+1, avg[pmap.index(i).hv()]);
+        hrms.SetBinContent(i+1, rms[pmap.index(i).hv()]);
+    }
+
+    return true;
+}
+
+bool extract_dac(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
+{
+    TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.BIAS_CONTROL_DAC.fits",
+                        day/10000, (day/100)%100, day%100, day);
+
+    // =====================================================================
+
+    // Now we read the temperatures of the sensors during the
+    // single pe run and average them
+    fits faux(naux.Data());
+    if (!faux)
+    {
+        cout << "Could not open " << naux << endl;
+        return false;
+    }
+
+    Double_t time;
+    UShort_t dac[416];
+
+    if (!faux.SetPtrAddress("Time",  &time))
+        return false;
+    if (!faux.SetPtrAddress("U", dac))
+        return false;
+
+    TArrayD avg(320);
+    TArrayD rms(320);
+    TArrayI cnt(320);
+
+    while (faux.GetNextRow())
+    {
+        if (time<tstart || time>tstop)
+            continue;
+
+        TArrayD med(320);
+
+        for (int i=0; i<320; i++)
+        {
+            med[i]  = dac[i];
+            avg[i] += dac[i];
+            rms[i] += dac[i]*dac[i];
+            cnt[i]++;
+        }
+
+        Double_t median;
+        MMath::MedianDev(320, med.GetArray(), median);
+        g.SetPoint(g.GetN(), time, median);
+    }
+
+    if (cnt[0]==0)
+        return 2;
+
+    for (int i=0; i<320; i++)
+    {
+        avg[i] /= cnt[i];
+        rms[i] /= cnt[i];
+        rms[i] -= avg[i]*avg[i];
+        rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
+    }
+
+    for (int i=0; i<1440; i++)
+    {
+        havg.SetBinContent(i+1, avg[pmap.index(i).hv()]);
+        hrms.SetBinContent(i+1, rms[pmap.index(i).hv()]);
+    }
+
+    return true;
+}
+
+bool extract_temp(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
+{
+    TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FSC_CONTROL_TEMPERATURE.fits",
+                        day/10000, (day/100)%100, day%100, day);
+
+    Interpolator2D interpol;
+
+    const auto sens = Interpolator2D::ReadGrid("operation/sensor-pos.txt");
+    if (sens.size()!=31)
+    {
+        cout << "Reading sensor-pos.txt failed: " << interpol.getInputGrid().size() << endl;
+        return false;
+    }
+
+    // =====================================================================
+
+    // Now we read the temperatures of the sensors during the
+    // single pe run and average them
+    fits faux(naux.Data());
+    if (!faux)
+    {
+        cout << "Could not open " << naux << endl;
+        return false;
+    }
+
+    Double_t time;
+    Float_t  temp[31];
+
+    if (!faux.SetPtrAddress("Time",  &time))
+        return false;
+    if (!faux.SetPtrAddress("T_sens", temp))
+        return false;
+
+    TArrayD avg(31);
+    TArrayD rms(31);
+    TArrayI cnt(31);
+
+    while (faux.GetNextRow())
+    {
+        if (time<tstart || time>tstop)
+            continue;
+
+        double mean  = 0;
+        int    count = 0;
+
+        for (int i=0; i<31; i++)
+        {
+            if (temp[i]==0)
+                continue;
+
+            double T = temp[i];
+
+            mean += T;
+            count++;
+
+            avg[i] += T;
+            rms[i] += T*T;
+            cnt[i]++;
+        }
+
+        g.SetPoint(g.GetN(), time, mean/count);
+
+    }
+
+    vector<double> zavg;
+    vector<double> zrms;
+
+    vector<Interpolator2D::vec> pos;
+    for (int i=0; i<31; i++)
+    {
+        if (cnt[i]>0)
+        {
+            avg[i] /= cnt[i];
+            rms[i] /= cnt[i];
+            rms[i] -= avg[i]*avg[i];
+            rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
+
+            pos.push_back(sens[i]);
+            zavg.push_back(avg[i]);
+            zrms.push_back(rms[i]);
+        }
+    }
+
+    // ================================================================
+
+    interpol.SetInputGrid(pos);
+
+    if (!interpol.ReadOutputGrid("operation/bias-positions.txt"))
+    {
+        cout << "Error setting output grid." << endl;
+        return false;
+    }
+    if (interpol.getOutputGrid().size()!=320)
+    {
+        cout << "Reading bias-positions.txt failed: " << interpol.getOutputGrid().size() << endl;
+        return false;
+    }
+
+    // ================================================================
+
+    const auto &input = interpol.getInputGrid();
+
+    TGraph sensors(input.size());
+
+    for (unsigned int i=0; i<input.size(); i++)
+        sensors.SetPoint(i, input[i].x, input[i].y);
+
+    const vector<double> rc_avg = interpol.Interpolate(zavg);
+    const vector<double> rc_rms = interpol.Interpolate(zrms);
+
+    // ================================================================
+
+    for (int i=0; i<1440; i++)
+    {
+        havg.SetBinContent(i+1, rc_avg[pmap.index(i).hv()]);
+        hrms.SetBinContent(i+1, rc_rms[pmap.index(i).hv()]);
+    }
+
+    return true;
+
+}
+
+void extract_aux(UInt_t day=20131114, UInt_t run=136, const char *output=0)
+{
+    if (!pmap.Read("operation/FACTmap111030.txt"))
+    {
+        cout << "FACTmap111030.txt not found." << endl;
+        return;
+    }
+
+    // =====================================================================
+
+    TString nraw = Form("/fact/raw/%4d/%02d/%02d/%8d_%03d.fits.fz",
+                        day/10000, (day/100)%100, day%100, day, run);
+
+
+    zfits fraw(nraw.Data());
+    if (!fraw)
+    {
+        cout << "Open raw file failed: " << nraw << endl;
+        return;
+    }
+
+    double tstart = fraw.GetUInt("TSTARTI") + fraw.GetFloat("TSTARTF");
+    double tstop  = fraw.GetUInt("TSTOPI")  + fraw.GetFloat("TSTOPF");
+
+    // =====================================================================
+
+    MGeomCamFACT geom;
+
+    TGraph gcur;
+    MHCamera hcur_m(geom);
+    MHCamera hcur_r(geom);
+    if (!extract_current(day, tstart, tstop, hcur_m, hcur_r, gcur))
+        return;
+
+    TGraph gtmp;
+    MHCamera htmp_m(geom);
+    MHCamera htmp_r(geom);
+    if (!extract_temp(day, tstart, tstop, htmp_m, htmp_r, gtmp))
+        return;
+
+    TGraph gdac;
+    MHCamera hdac_m(geom);
+    MHCamera hdac_r(geom);
+    if (!extract_dac(day, tstart, tstop, hdac_m, hdac_r, gdac))
+        return;
+
+    // ================================================================
+
+    TCanvas *c = new TCanvas;
+    c->Divide(3, 3);
+
+    c->cd(1);
+    hcur_m.SetName("CurrentAvg");
+    hcur_m.SetAllUsed();
+    hcur_m.DrawCopy();
+    c->cd(4);
+    hcur_r.SetName("CurrentRms");
+    hcur_r.SetAllUsed();
+    hcur_r.DrawCopy();
+    c->cd(7);
+    gcur.SetName("CurrentTime");
+    gcur.SetMarkerStyle(kFullDotMedium);
+    gcur.SetMinimum(0);
+    gcur.DrawClone("AP");
+
+    c->cd(2);
+    htmp_m.SetName("TempAvg");
+    htmp_m.SetAllUsed();
+    htmp_m.DrawCopy();
+    c->cd(5);
+    htmp_r.SetName("TempRms");
+    htmp_r.SetAllUsed();
+    htmp_r.DrawCopy();
+    c->cd(8);
+    gtmp.SetName("TempTime");
+    gtmp.SetMarkerStyle(kFullDotMedium);
+    gtmp.SetMinimum(0);
+    gtmp.DrawClone("AP");
+
+    c->cd(3);
+    hdac_m.SetName("DacAvg");
+    hdac_m.SetAllUsed();
+    hdac_m.DrawCopy();
+    c->cd(6);
+    hdac_r.SetName("DacRms");
+    hdac_r.SetAllUsed();
+    hdac_r.DrawCopy();
+    c->cd(9);
+    gdac.SetName("DacTime");
+    gdac.SetMarkerStyle(kFullDotMedium);
+    gdac.SetMinimum(0);
+    gdac.DrawClone("AP");
+
+    if (output)
+        c->SaveAs(output);
+}