Changeset 17370 for trunk


Ignore:
Timestamp:
11/25/13 18:09:42 (11 years ago)
Author:
tbretz
Message:
Renamed and added more auxiliary infomrations.
File:
1 copied

Legend:

Unmodified
Added
Removed

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

    r17326 r17370  
    99#include "Interpolator2D.h"
    1010
     11#include "MMath.h"
    1112#include "MHCamera.h"
    1213#include "MGeomCamFACT.h"
    1314
    14 int extract_temp(UInt_t day=20130902, UInt_t run=146, const char *output=0)
     15PixelMap pmap;
     16
     17bool extract_current(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
    1518{
    16     TString nraw = Form("/fact/raw/%4d/%02d/%02d/%8d_%03d.fits.fz",
    17                         day/10000, (day/100)%100, day%100, day, run);
    18 
    19     TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FSC_CONTROL_TEMPERATURE.fits",
     19    TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FEEDBACK_CALIBRATED_CURRENTS.fits",
    2020                        day/10000, (day/100)%100, day%100, day);
    21 
    22     PixelMap pmap;
    23     if (!pmap.Read("FACTmap111030.txt"))
    24     {
    25         cout << "FACTmap111030.txt not found." << endl;
    26         return 1;
    27     }
    28 
    29     Interpolator2D interpol;
    30 
    31     const auto sens = Interpolator2D::ReadGrid("operation/sensor-pos.txt");
    32     if (sens.size()!=31)
    33     {
    34         cout << "Reading sensor-pos.txt failed: " << interpol.getInputGrid().size() << endl;
    35         return 2;
    36     }
    37 
    38     // =====================================================================
    39 
    40     zfits fraw(nraw.Data());
    41     if (!fraw)
    42     {
    43         cout << "Open raw file failed: " << nraw << endl;
    44         return 1;
    45     }
    46 
    47     double tstart = fraw.GetUInt("TSTARTI") + fraw.GetFloat("TSTARTF");
    48     double tstop  = fraw.GetUInt("TSTOPI")  + fraw.GetFloat("TSTOPF");
    4921
    5022    // =====================================================================
    … …  
    5628    {
    5729        cout << "Could not open " << naux << endl;
    58         return 1;
     30        return false;
     31    }
     32
     33    Double_t time;
     34    Float_t  I[416];
     35
     36    if (!faux.SetPtrAddress("Time",  &time))
     37        return false;
     38    if (!faux.SetPtrAddress("I", I))
     39        return false;
     40
     41    TArrayD avg(320);
     42    TArrayD rms(320);
     43    TArrayI cnt(320);
     44
     45    while (faux.GetNextRow())
     46    {
     47        if (time<tstart || time>tstop)
     48            continue;
     49
     50        TArrayD med(320);
     51
     52        double mean = 0;
     53
     54        int j=0;
     55        for (int i=0; i<320; i++)
     56        {
     57            if (i!=66 && i!=191 && i!=193 && i!=17 && i!=206)
     58                med[j++]  = I[i];
     59
     60            avg[i] += I[i];
     61            rms[i] += I[i]*I[i];
     62            cnt[i]++;
     63        }
     64
     65        Double_t median;
     66        MMath::MedianDev(315, med.GetArray(), median);
     67        g.SetPoint(g.GetN(), time, median);
     68    }
     69
     70    if (cnt[0]==0)
     71        return 2;
     72
     73    for (int i=0; i<320; i++)
     74    {
     75        avg[i] /= cnt[i];
     76        rms[i] /= cnt[i];
     77        rms[i] -= avg[i]*avg[i];
     78        rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
     79    }
     80
     81    double mn = 0;
     82    for (int i=0; i<1440; i++)
     83    {
     84        mn += avg[pmap.index(i).hv()];
     85        havg.SetBinContent(i+1, avg[pmap.index(i).hv()]);
     86        hrms.SetBinContent(i+1, rms[pmap.index(i).hv()]);
     87    }
     88
     89    return true;
     90}
     91
     92bool extract_dac(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
     93{
     94    TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.BIAS_CONTROL_DAC.fits",
     95                        day/10000, (day/100)%100, day%100, day);
     96
     97    // =====================================================================
     98
     99    // Now we read the temperatures of the sensors during the
     100    // single pe run and average them
     101    fits faux(naux.Data());
     102    if (!faux)
     103    {
     104        cout << "Could not open " << naux << endl;
     105        return false;
     106    }
     107
     108    Double_t time;
     109    UShort_t dac[416];
     110
     111    if (!faux.SetPtrAddress("Time",  &time))
     112        return false;
     113    if (!faux.SetPtrAddress("U", dac))
     114        return false;
     115
     116    TArrayD avg(320);
     117    TArrayD rms(320);
     118    TArrayI cnt(320);
     119
     120    while (faux.GetNextRow())
     121    {
     122        if (time<tstart || time>tstop)
     123            continue;
     124
     125        TArrayD med(320);
     126
     127        for (int i=0; i<320; i++)
     128        {
     129            med[i]  = dac[i];
     130            avg[i] += dac[i];
     131            rms[i] += dac[i]*dac[i];
     132            cnt[i]++;
     133        }
     134
     135        Double_t median;
     136        MMath::MedianDev(320, med.GetArray(), median);
     137        g.SetPoint(g.GetN(), time, median);
     138    }
     139
     140    if (cnt[0]==0)
     141        return 2;
     142
     143    for (int i=0; i<320; i++)
     144    {
     145        avg[i] /= cnt[i];
     146        rms[i] /= cnt[i];
     147        rms[i] -= avg[i]*avg[i];
     148        rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
     149    }
     150
     151    for (int i=0; i<1440; i++)
     152    {
     153        havg.SetBinContent(i+1, avg[pmap.index(i).hv()]);
     154        hrms.SetBinContent(i+1, rms[pmap.index(i).hv()]);
     155    }
     156
     157    return true;
     158}
     159
     160bool extract_temp(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
     161{
     162    TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FSC_CONTROL_TEMPERATURE.fits",
     163                        day/10000, (day/100)%100, day%100, day);
     164
     165    Interpolator2D interpol;
     166
     167    const auto sens = Interpolator2D::ReadGrid("operation/sensor-pos.txt");
     168    if (sens.size()!=31)
     169    {
     170        cout << "Reading sensor-pos.txt failed: " << interpol.getInputGrid().size() << endl;
     171        return false;
     172    }
     173
     174    // =====================================================================
     175
     176    // Now we read the temperatures of the sensors during the
     177    // single pe run and average them
     178    fits faux(naux.Data());
     179    if (!faux)
     180    {
     181        cout << "Could not open " << naux << endl;
     182        return false;
    59183    }
    60184
    … …  
    63187
    64188    if (!faux.SetPtrAddress("Time",  &time))
    65         return -1;
     189        return false;
    66190    if (!faux.SetPtrAddress("T_sens", temp))
    67         return -1;
     191        return false;
    68192
    69193    TArrayD avg(31);
    … …  
    75199        if (time<tstart || time>tstop)
    76200            continue;
     201
     202        double mean  = 0;
     203        int    count = 0;
    77204
    78205        for (int i=0; i<31; i++)
    … …  
    83210            double T = temp[i];
    84211
     212            mean += T;
     213            count++;
     214
    85215            avg[i] += T;
    86216            rms[i] += T*T;
    87217            cnt[i]++;
    88218        }
    89     }
    90 
    91     vector<double> z;
     219
     220        g.SetPoint(g.GetN(), time, mean/count);
     221
     222    }
     223
     224    vector<double> zavg;
     225    vector<double> zrms;
    92226
    93227    vector<Interpolator2D::vec> pos;
    … …  
    98232            avg[i] /= cnt[i];
    99233            rms[i] /= cnt[i];
    100             rms[i] = sqrt(rms[i]-avg[i]*avg[i]);
     234            rms[i] -= avg[i]*avg[i];
     235            rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
    101236
    102237            pos.push_back(sens[i]);
    103             z.push_back(avg[i]);
     238            zavg.push_back(avg[i]);
     239            zrms.push_back(rms[i]);
    104240        }
    105241    }
    … …  
    112248    {
    113249        cout << "Error setting output grid." << endl;
    114         return 3;
     250        return false;
    115251    }
    116252    if (interpol.getOutputGrid().size()!=320)
    117253    {
    118254        cout << "Reading bias-positions.txt failed: " << interpol.getOutputGrid().size() << endl;
    119         return 3;
     255        return false;
    120256    }
    121257
    … …  
    129265        sensors.SetPoint(i, input[i].x, input[i].y);
    130266
    131     const vector<double> rc = interpol.Interpolate(z);
     267    const vector<double> rc_avg = interpol.Interpolate(zavg);
     268    const vector<double> rc_rms = interpol.Interpolate(zrms);
    132269
    133270    // ================================================================
    134271
     272    for (int i=0; i<1440; i++)
     273    {
     274        havg.SetBinContent(i+1, rc_avg[pmap.index(i).hv()]);
     275        hrms.SetBinContent(i+1, rc_rms[pmap.index(i).hv()]);
     276    }
     277
     278    return true;
     279
     280}
     281
     282void extract_aux(UInt_t day=20131114, UInt_t run=136, const char *output=0)
     283{
     284    if (!pmap.Read("operation/FACTmap111030.txt"))
     285    {
     286        cout << "FACTmap111030.txt not found." << endl;
     287        return;
     288    }
     289
     290    // =====================================================================
     291
     292    TString nraw = Form("/fact/raw/%4d/%02d/%02d/%8d_%03d.fits.fz",
     293                        day/10000, (day/100)%100, day%100, day, run);
     294
     295
     296    zfits fraw(nraw.Data());
     297    if (!fraw)
     298    {
     299        cout << "Open raw file failed: " << nraw << endl;
     300        return;
     301    }
     302
     303    double tstart = fraw.GetUInt("TSTARTI") + fraw.GetFloat("TSTARTF");
     304    double tstop  = fraw.GetUInt("TSTOPI")  + fraw.GetFloat("TSTOPF");
     305
     306    // =====================================================================
     307
     308    MGeomCamFACT geom;
     309
     310    TGraph gcur;
     311    MHCamera hcur_m(geom);
     312    MHCamera hcur_r(geom);
     313    if (!extract_current(day, tstart, tstop, hcur_m, hcur_r, gcur))
     314        return;
     315
     316    TGraph gtmp;
     317    MHCamera htmp_m(geom);
     318    MHCamera htmp_r(geom);
     319    if (!extract_temp(day, tstart, tstop, htmp_m, htmp_r, gtmp))
     320        return;
     321
     322    TGraph gdac;
     323    MHCamera hdac_m(geom);
     324    MHCamera hdac_r(geom);
     325    if (!extract_dac(day, tstart, tstop, hdac_m, hdac_r, gdac))
     326        return;
     327
     328    // ================================================================
     329
    135330    TCanvas *c = new TCanvas;
    136     c->SetName("Temp");
    137 
    138     MGeomCamFACT geom;
    139 
    140     MHCamera h1(geom);
    141     h1.SetName("Temp");
    142     for (int i=0; i<1440; i++)
    143         h1.SetBinContent(i+1, rc[pmap.index(i).hv()]);
    144 
    145     h1.SetAllUsed();
    146     h1.DrawCopy();
    147 
    148     sensors.SetName("Sensors");
    149     sensors.SetMarkerColor(kWhite);
    150     sensors.DrawClone("*");
     331    c->Divide(3, 3);
     332
     333    c->cd(1);
     334    hcur_m.SetName("CurrentAvg");
     335    hcur_m.SetAllUsed();
     336    hcur_m.DrawCopy();
     337    c->cd(4);
     338    hcur_r.SetName("CurrentRms");
     339    hcur_r.SetAllUsed();
     340    hcur_r.DrawCopy();
     341    c->cd(7);
     342    gcur.SetName("CurrentTime");
     343    gcur.SetMarkerStyle(kFullDotMedium);
     344    gcur.SetMinimum(0);
     345    gcur.DrawClone("AP");
     346
     347    c->cd(2);
     348    htmp_m.SetName("TempAvg");
     349    htmp_m.SetAllUsed();
     350    htmp_m.DrawCopy();
     351    c->cd(5);
     352    htmp_r.SetName("TempRms");
     353    htmp_r.SetAllUsed();
     354    htmp_r.DrawCopy();
     355    c->cd(8);
     356    gtmp.SetName("TempTime");
     357    gtmp.SetMarkerStyle(kFullDotMedium);
     358    gtmp.SetMinimum(0);
     359    gtmp.DrawClone("AP");
     360
     361    c->cd(3);
     362    hdac_m.SetName("DacAvg");
     363    hdac_m.SetAllUsed();
     364    hdac_m.DrawCopy();
     365    c->cd(6);
     366    hdac_r.SetName("DacRms");
     367    hdac_r.SetAllUsed();
     368    hdac_r.DrawCopy();
     369    c->cd(9);
     370    gdac.SetName("DacTime");
     371    gdac.SetMarkerStyle(kFullDotMedium);
     372    gdac.SetMinimum(0);
     373    gdac.DrawClone("AP");
    151374
    152375    if (output)
    153376        c->SaveAs(output);
    154 
    155     return 0;
    156377}
Note: See TracChangeset for help on using the changeset viewer.