Changeset 17370 for trunk/Mars/fact/analysis/gain/extract_aux.C

Timestamp:

11/25/13 18:09:42 (11 years ago)

Author:

tbretz

Message:

Renamed and added more auxiliary infomrations.

File:

• trunk/Mars/fact/analysis/gain/extract_aux.C (copied) (copied from trunk/Mars/fact/analysis/gain/extract_temp.C ) (8 diffs)

trunk/Mars/fact/analysis/gain/extract_aux.C

@@ -9 +9 @@
 #include "Interpolator2D.h"
 
+#include "MMath.h"
 #include "MHCamera.h"
 #include "MGeomCamFACT.h"
 
-int extract_temp(UInt_t day=20130902, UInt_t run=146, const char *output=0)
+PixelMap pmap;
+
+bool extract_current(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
 {
-    TString nraw = Form("/fact/raw/%4d/%02d/%02d/%8d_%03d.fits.fz",
-                        day/10000, (day/100)%100, day%100, day, run);
-
-    TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FSC_CONTROL_TEMPERATURE.fits",
+    TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FEEDBACK_CALIBRATED_CURRENTS.fits",
                         day/10000, (day/100)%100, day%100, day);
-
-    PixelMap pmap;
-    if (!pmap.Read("FACTmap111030.txt"))
-    {
-        cout << "FACTmap111030.txt not found." << endl;
-        return 1;
-    }
-
-    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 2;
-    }
-
-    // =====================================================================
-
-    zfits fraw(nraw.Data());
-    if (!fraw)
-    {
-        cout << "Open raw file failed: " << nraw << endl;
-        return 1;
-    }
-
-    double tstart = fraw.GetUInt("TSTARTI") + fraw.GetFloat("TSTARTF");
-    double tstop  = fraw.GetUInt("TSTOPI")  + fraw.GetFloat("TSTOPF");
 
     // =====================================================================
@@ -56 +28 @@
     {
         cout << "Could not open " << naux << endl;
-        return 1;
+        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;
     }
 
@@ -63 +187 @@
 
     if (!faux.SetPtrAddress("Time",  &time))
-        return -1;
+        return false;
     if (!faux.SetPtrAddress("T_sens", temp))
-        return -1;
+        return false;
 
     TArrayD avg(31);
@@ -75 +199 @@
         if (time<tstart || time>tstop)
             continue;
+
+        double mean  = 0;
+        int    count = 0;
 
         for (int i=0; i<31; i++)
@@ -83 +210 @@
             double T = temp[i];
 
+            mean += T;
+            count++;
+
             avg[i] += T;
             rms[i] += T*T;
             cnt[i]++;
         }
-    }
-
-    vector<double> z;
+
+        g.SetPoint(g.GetN(), time, mean/count);
+
+    }
+
+    vector<double> zavg;
+    vector<double> zrms;
 
     vector<Interpolator2D::vec> pos;
@@ -98 +232 @@
             avg[i] /= cnt[i];
             rms[i] /= cnt[i];
-            rms[i] = sqrt(rms[i]-avg[i]*avg[i]);
+            rms[i] -= avg[i]*avg[i];
+            rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
 
             pos.push_back(sens[i]);
-            z.push_back(avg[i]);
+            zavg.push_back(avg[i]);
+            zrms.push_back(rms[i]);
         }
     }
@@ -112 +248 @@
     {
         cout << "Error setting output grid." << endl;
-        return 3;
+        return false;
     }
     if (interpol.getOutputGrid().size()!=320)
     {
         cout << "Reading bias-positions.txt failed: " << interpol.getOutputGrid().size() << endl;
-        return 3;
+        return false;
     }
 
@@ -129 +265 @@
         sensors.SetPoint(i, input[i].x, input[i].y);
 
-    const vector<double> rc = interpol.Interpolate(z);
+    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->SetName("Temp");
-
-    MGeomCamFACT geom;
-
-    MHCamera h1(geom);
-    h1.SetName("Temp");
-    for (int i=0; i<1440; i++)
-        h1.SetBinContent(i+1, rc[pmap.index(i).hv()]);
-
-    h1.SetAllUsed();
-    h1.DrawCopy();
-
-    sensors.SetName("Sensors");
-    sensors.SetMarkerColor(kWhite);
-    sensors.DrawClone("*");
+    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);
-
-    return 0;
 }