#include <algorithm>

void threshold(const char *fname, double beg=0, double end=100000)
{
    fits file(fname);

    //file.PrintColumns();
    //file.PrintKeys();

    Double_t time;
    UShort_t adc;
    UInt_t   dead, run;
    float    temp;
    UInt_t   qos;
    UInt_t   counter_ch[8];
    UInt_t   counter_sum[2];
    UShort_t dac_ch[8];
    UShort_t dac_sum[2];
    float    dt_sec[2];
    file.SetPtrAddress("QoS",         &qos);
    file.SetPtrAddress("Time",        &time);
    file.SetPtrAddress("counter_ch",   counter_ch);
    file.SetPtrAddress("counter_sum",  counter_sum);
    file.SetPtrAddress("dac_ch",       dac_ch);
    file.SetPtrAddress("dac_sum",      dac_sum);
    file.SetPtrAddress("dt_sec",       dt_sec);

    UInt_t offset = file.GetUInt("MJDREF");
    if (beg < 30000)
        beg+=offset;
    if (end < 30000)
        end+=offset;

    double min_dac_ch   = 65535;
    double avg_dac_ch   = 0;
    double max_dac_ch   = 0;

    double dac_ch_min   = 65535;
    double dac_ch_med   = 0;
    double dac_ch_max   = 0;

    double min_dac_sum  = 65535;
    double max_dac_sum  = 0;

    double dac_sum_min  = 65535;
    double dac_sum_max  = 0;

    double rate_ch[8] = {0,0,0,0,0,0,0,0};
    double rate_sum[2] = {0,0};

    int count = 0;

    double diff = -1;

    double sum_dt = 0;


    while (file.GetNextRow())
    {
        time += offset;

        if (time>end)
            break;

        dac_ch_min  = TMath::MinElement(8, dac_ch);
        dac_ch_med  = TMath::Median(    8, dac_ch);
        dac_ch_max  = TMath::MaxElement(8, dac_ch);
        dac_sum_min = TMath::Min(dac_sum[0], dac_sum[1]);
        dac_sum_max = TMath::Min(dac_sum[0], dac_sum[1]);

        if (time<beg)
        {
            diff = beg-time;
            continue;
        }

        avg_dac_ch += dac_ch_med;
        count++;

        min_dac_ch = TMath::Min(min_dac_ch, dac_ch_min);
        max_dac_ch = TMath::Max(max_dac_ch, dac_ch_max);

        min_dac_sum = TMath::Min(min_dac_sum, dac_sum_min);
        max_dac_sum = TMath::Max(max_dac_sum, dac_sum_max);

        for (int i=0;i<8; i++)
            rate_ch[i] += counter_ch[i];
        rate_sum[0] += counter_sum[0];
        rate_sum[1] += counter_sum[1];

        sum_dt += dt_sec[0];
    }

    if (count==0)
    {
        if (diff<5./24/3600)
            return;

        cout << "result " << dac_ch_min << ", " << dac_ch_med << ", " << dac_ch_max << ", " << dac_sum_min << ", " << dac_sum_max;
        cout << ", NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL" << endl;
        return;
    }

    avg_dac_ch /= count;


    cout << "result " << min_dac_ch << ", " << avg_dac_ch << ", " << max_dac_ch << ", " << min_dac_sum << ", " << max_dac_sum;
    for (int i=0;i<8; i++)
        cout << ", " << rate_ch[i] / sum_dt;
    cout << ", " << rate_sum[0] / sum_dt;
    cout << ", " << rate_sum[1] / sum_dt;
    cout << endl;
}