int StateMachineFeedback::HandleBiasCurrent ( const EventImp &  evt )

inlineprivate

Definition at line 457 of file feedback.cc.

References __attribute__, begin, PixelMapEntry::count(), data, Debug, EventImp::GetSize(), EventImp::GetTime(), PixelMapEntry::group(), PixelMap::hv(), i, Time::IsValid(), Feedback::State::kCalibrated, Feedback::State::kCalibrating, Feedback::State::kConnected, Feedback::State::kInProgress, BIAS::kNumChannels, Feedback::State::kOnStandby, BIAS::State::kRamping, BIAS::State::kVoltageOn, Feedback::State::kWaitingForData, EventImp::Ptr(), Dim::SendCommandNB(), DimClient::sendCommandNB(), DimDescribedService::setData(), DimService::setQuality(), DimState::state(), Time::UnixTime(), and DimDescribedService::Update().

Referenced by StateMachineFeedback().

    {
        if (!CheckEventSize(evt.GetSize(), "HandleBiasCurrent", BIAS::kNumChannels*sizeof(uint16_t)))
            return Feedback::State::kConnected;

        if (GetCurrentState()<Feedback::State::kCalibrating)
            return GetCurrentState();

        // ------------------------------- HandleCalibration -----------------------------------
        if (GetCurrentState()==Feedback::State::kCalibrating)
            return HandleCalibration(evt);

        // ---------------------- Calibrated, WaitingForData, InProgress -----------------------

        // We are waiting but no valid temperature yet, go on waiting
        if (GetCurrentState()==Feedback::State::kWaitingForData &&
            (!fTimeTemp.IsValid() || Time()-fTimeTemp>boost::posix_time::minutes(5)))
            return GetCurrentState();

        // We are waiting but biasctrl is still in ramping (this might
        // be the case if the feedback was started with a new overvoltage
        // while the last ramping command was still in progress)
        if (GetCurrentState()==Feedback::State::kWaitingForData &&
            fDimBias.state()==BIAS::State::kRamping)
            return GetCurrentState();

        // We are already in progress but no valid temperature update anymore
        if (GetCurrentState()>=Feedback::State::kInProgress &&
            (!fTimeTemp.IsValid() || Time()-fTimeTemp>boost::posix_time::minutes(5)))
        {
            Warn("Current control in progress, but last received temperature older than 5min... switching voltage off.");
            Dim::SendCommandNB("BIAS_CONTROL/SET_ZERO_VOLTAGE");
            return Feedback::State::kCalibrated;
        }

        // ---------------------- Calibrated, WaitingForData, InProgress -----------------------

        const int Navg = fDimBias.state()!=BIAS::State::kVoltageOn ? 1 : 3;

        const vector<float> &Imes = AverageCurrents(evt.Ptr<int16_t>(), Navg).first;
        if (Imes.size()==0)
            return GetCurrentState();

        fCurrentsAvg.assign(BIAS::kNumChannels, 0);
        fCurrentsRms.assign(BIAS::kNumChannels, 0);
        fCursorCur = 0;

        // -------------------------------------------------------------------------------------
        // Inner patches to be blocked (operated below the operation voltage) in moon mode

        static const array<int, 14> inner0 =
        {{
             62,  63, 130, 131, 132, 133, 134,
            135, 222, 223, 292, 293, 294, 295,
        }};

        static const array<int, 23> inner1 =
        {{
             58,  59,  60,  61, 129, 138, 139, 140, 141, 142, 143, 218,
            219, 220, 221, 290, 291, 298, 299, 300, 301, 302, 303,
        }};

        static const array<int, 43> inner2 =
        {{
             42,  43,  44,  45,  55,  56,  57,  70,  71,  78,  79,
             96,  97,  98,  99, 102, 103, 128, 136, 137, 159, 202,
            203, 204, 205, 214, 216, 217, 228, 230, 231, 256, 257,
            258, 259, 262, 263, 288, 289, 296, 297, 310, 318
        }};

        // -------------------------------------------------------------------------------------

        // Nominal overvoltage (w.r.t. the bias setup values)
        const double voltageoffset = GetCurrentState()<Feedback::State::kWaitingForData ? 0 : fUserOffset;

        double avg[2] = {   0,   0 };
        double min[2] = {  90,  90 };
        double max[2] = { -90, -90 };
        int    num[3] = {   0,   0,   0 };

        vector<double> med[3];
        med[0].resize(BIAS::kNumChannels);
        med[1].resize(BIAS::kNumChannels);
        med[2].resize(BIAS::kNumChannels);

        struct dim_data
        {
            float I[BIAS::kNumChannels];
            float Iavg;
            float Irms;
            float Imed;
            float Idev;
            uint32_t N;
            float Tdiff;
            float Uov[BIAS::kNumChannels];
            float Unom;
            float dUtemp;

            dim_data() { memset(this, 0, sizeof(dim_data)); }
        } __attribute__((__packed__));

        int Ndev[3] = { 0, 0, 0 };

        dim_data data;

        data.Unom   = voltageoffset;
        data.dUtemp = fTempOffsetAvg;

        vector<float> vec(BIAS::kNumChannels);

        // ================================= old =======================
        // Pixel  583: 5 31 == 191 (5)  C2 B3 P3
        // Pixel  830: 2  2 ==  66 (4)  C0 B8 P1
        // Pixel 1401: 6  1 == 193 (5)  C2 B4 P0

        double UdrpAvg = 0;
        double UdrpRms = 0;

        for (int i=0; i<320/*BIAS::kNumChannels*/; i++)
        {
            const PixelMapEntry &hv = fMap.hv(i);
            if (!hv)
                continue;

            // Check if this is a blocked channel
            // 272 is the one with the shortcut
            const bool blocked =
                (fMoonMode>0 && std::find(inner0.begin(), inner0.end(), i)!=inner0.end()) ||
                (fMoonMode>1 && std::find(inner1.begin(), inner1.end(), i)!=inner1.end()) ||
                (fMoonMode>2 && std::find(inner2.begin(), inner2.end(), i)!=inner2.end()) ||
                i==272;

            // Number of G-APDs in this patch
            const int N = hv.count();

            // Average measured ADC value for this channel
            // FIXME: This is a workaround for the problem with the
            // readout of bias voltage channel 263
            const double adc = Imes[i]/* * (5e-3/4096)*/; // [A]

            // Current through ~100 Ohm measurement resistor
            //const double I8 = (adc-fCalibDeltaI[i])*fCalibR8[i]/100;
            const double I8 = adc-fCalibDeltaI[i];

            // Current through calibration resistors (R9)
            // This is uncalibrated, but since the corresponding calibrated
            // value I8 is subtracted, the difference should yield a correct value
            const double I9 = fBiasDac[i] * (1e-3/4096);//U9/R9;   [A]

            // Current in R4/R5 branch
            //const double Iout = I8 - I9;//I8>I9 ? I8 - I9 : 0;
            const double Iout = I8 - I9*100/fCalibR8[i];//I8>I9 ? I8 - I9 : 0;

            // Applied voltage at calibration resistors, according to biasctrl
            const double U9 = fBiasVolt[i];

            //          new    I8 - I9*100/fCalibR8       100
            // change = --- = ---------------------- =  --------  = 0.8
            //          old    I8*fCalibR8/100 - I9     fCalibR8

            // Serial resistors (one 1kOhm at the output of the bias crate, one 1kOhm in the camera)
            const double R4 = 2000;

            // Serial resistor of the individual G-APDs plus 50 Ohm termination
            double R5 = 3900./N + 50;

            // This is assuming that the broken pixels have a 390 Ohm instead of 3900 Ohm serial resistor
            if (i==66 || i==193)               // Pixel 830(66) / Pixel 583(191)
                R5 = 1./((N-1)/3900.+1/1000.);
            if (i==191)                        // Pixel 1399(193)
                R5 = 1./((N-1)/3900.+1/390.);
            if (i==17 || i==206)               // dead pixel 923(80) / dead pixel 424(927)
                R5 = 3900./(N-1);              // cannot identify third dead pixel in light-pulser data

            // The measurement resistor
            const double R8 = 0;

            // Total resistance of branch with diodes (R4+R5)
            // Assuming that the voltage output of the OpAMP is linear
            // with the DAC setting and not the voltage at R9, the
            // additional voltage drop at R8 must be taken into account
            const double R = R4 + R5 + R8;

            // For the patches with a broken resistor - ignoring the G-APD resistance -
            // we get:
            //
            // I[R=3900] =  Iout *      1/(10+(N-1))  = Iout        /(N+9)
            // I[R= 390] =  Iout * (1 - 1/(10+(N-1))) = Iout * (N+8)/(N+9)
            //
            // I[R=390] / I[R=3900] = N+8
            //
            // Udrp = Iout*3900/(N+9) + Iout*1000 + Iout*1000 = Iout * R

            // Voltage drop in R4/R5 branch (for the G-APDs with correct resistor)
            // The voltage drop should not be <0, otherwise an unphysical value
            // would be amplified when Uset is calculated.
            const double Udrp = Iout<0 ? 0 : R*Iout;

            // Nominal operation voltage with correction for temperature dependence
            const double Uop = fVoltGapd[i] + fVoltOffset[i] + fTempOffset[i]
                + (blocked ? -5 : 0);

            // Current overvoltage (at a G-APD with the correct 3900 Ohm resistor)
            // expressed w.r.t. to the operation voltage
            const double Uov = (U9-Udrp)-Uop>-1.4 ? (U9-Udrp)-Uop : -1.4;

            // The current through one G-APD is the sum divided by the number of G-APDs
            // (assuming identical serial resistors)
            double Iapd = Iout/N;

            // Rtot = Uapd/Iout
            // Ich  = Uapd/Rch = (Rtot*Iout) / Rch = Rtot/Rch * Iout
            //
            // Rtot = 3900/N
            // Rch  = 3900
            //
            // Rtot = 1./((N-1)/3900 + 1/X)       X=390 or X=1000
            // Rch  = 3900
            //
            // Rtot/Rch =   1/((N-1)/3900 + 1/X)/3900
            // Rtot/Rch =   1/( [ X*(N-1) + 3900 ] / [ 3900 * X ])/3900
            // Rtot/Rch =   X/( [ X*(N-1)/3900 + 1 ] )/3900
            // Rtot/Rch =   X/( [ X*(N-1) + 3900 ] )
            // Rtot/Rch =   1/( [ (N-1) + 3900/X ] )
            //
            // Rtot/Rch[390Ohm]  =  1/( [ N + 9.0 ] )
            // Rtot/Rch[1000Ohm] =  1/( [ N + 2.9 ] )
            //
            // In this and the previosu case we neglect the resistance of the G-APDs, but we can make an
            // assumption: The differential resistance depends more on the NSB than on the PDE,
            // thus it is at least comparable for all G-APDs in the patch. In addition, although the
            // G-APD with the 390Ohm serial resistor has the wrong voltage applied, this does not
            // significantly influences the ohmic resistor or the G-APD because the differential
            // resistor is large enough that the increase of the overvoltage does not dramatically
            // increase the current flow as compared to the total current flow.
            if (i==66 || i==193)           // Iout/13 15.8   / Iout/14  16.8
                Iapd = Iout/(N+2.9);
            if (i==191)                    // Iout/7.9  38.3
                Iapd = Iout/(N+9);
            if (i==17 || i==206)
                Iapd = Iout/(N-1);

            // The differential resistance of the G-APD, i.e. the dependence of the
            // current above the breakdown voltage, is given by
            //const double Rapd = Uov/Iapd;
            // This allows us to estimate the current Iov at the overvoltage we want to apply
            //const double Iov = overvoltage/Rapd;

            // Estimate set point for over-voltage (voltage drop at the target point)
            // This estimation is based on the linear increase of the
            // gain with voltage and the increase of the crosstalk with
            // voltage, as measured with the overvoltage-tests (OVTEST)
            /*
             Uov+0.44<0.022 ?
                Ubd + overvoltage + Udrp*exp(0.6*(overvoltage-Uov))*pow((overvoltage+0.44), 0.6) :
                Ubd + overvoltage + Udrp*exp(0.6*(overvoltage-Uov))*pow((overvoltage+0.44)/(Uov+0.44), 0.6);
             */
            const double Uset =
                Uov+1.4<0.022 ?
                Uop + voltageoffset + Udrp*exp(0.6*(voltageoffset-Uov))*pow((voltageoffset+1.4),           0.6) :
                Uop + voltageoffset + Udrp*exp(0.6*(voltageoffset-Uov))*pow((voltageoffset+1.4)/(Uov+1.4), 0.6);

            if (fabs(voltageoffset-Uov)>0.033)
                Ndev[0]++;
            if (fabs(voltageoffset-Uov)>0.022)
                Ndev[1]++;
            if (fabs(voltageoffset-Uov)>0.011)
                Ndev[2]++;

            // Voltage set point
            vec[i] = Uset;

            const double iapd = Iapd*1e6; // A --> uA

            data.I[i]   = iapd;
            data.Uov[i] = Uov;

            if (!blocked)
            {
                const int g = hv.group();

                med[g][num[g]] = Uov;
                avg[g] += Uov;
                num[g]++;

                if (Uov<min[g])
                    min[g] = Uov;
                if (Uov>max[g])
                    max[g] = Uov;

                data.Iavg += iapd;
                data.Irms += iapd*iapd;

                med[2][num[2]++] = iapd;

                UdrpAvg += Udrp;
                UdrpRms += Udrp*Udrp;
            }
        }


        // ---------------------------- Calculate statistics ----------------------------------

        // average and rms
        data.Iavg /= num[2];
        data.Irms /= num[2];
        data.Irms -= data.Iavg*data.Iavg;

        data.N = num[2];
        data.Irms = data.Irms<0 ? 0: sqrt(data.Irms);

        // median
        sort(med[2].data(), med[2].data()+num[2]);

        data.Imed = num[2]%2 ? med[2][num[2]/2] : (med[2][num[2]/2-1]+med[2][num[2]/2])/2;

        // deviation
        for (int i=0; i<num[2]; i++)
            med[2][i] = fabs(med[2][i]-data.Imed);

        sort(med[2].data(), med[2].data()+num[2]);

        data.Idev = med[2][uint32_t(0.682689477208650697*num[2])];

        // time difference to calibration
        data.Tdiff = evt.GetTime().UnixTime()-fTimeCalib.UnixTime();

        // Average overvoltage
        const double Uov = (avg[0]+avg[1])/(num[0]+num[1]);

        // ------------------------------- Update voltages ------------------------------------

        int newstate = GetCurrentState();

        if (GetCurrentState()!=Feedback::State::kCalibrated) // WaitingForData, OnStandby, InProgress, kWarning, kCritical
        {
            if (fDimBias.state()!=BIAS::State::kRamping)
            {
                newstate = CheckLimits(data.I);

                // standby and change reduction level of voltage
                if (newstate==Feedback::State::kOnStandby)
                {
                    // Calculate average applied overvoltage and estimate an offset
                    // to reach fAbsoluteMedianCurrentLimit
                    float fAbsoluteMedianCurrentLimit = 85;
                    const double deltaU = (Uov+1.4)*(1-pow(fAbsoluteMedianCurrentLimit/data.Imed, 1./1.7));

                    if (fVoltageReduction+deltaU<0.033)
                        fVoltageReduction = 0;
                    else
                    {
                        fVoltageReduction += deltaU;

                        for (int i=0; i<320; i++)
                            vec[i] -= fVoltageReduction;
                    }
                }

                // FIXME: What if the brightest pixel gets too bright???
                // FIXME: What if fVolatgeReduction > U1.4V?

                // set voltage in 262 -> current in 262/263
                vec[263] = vec[262]-fVoltGapd[262]+fVoltGapd[263];

                // Do not ramp the channel with a shortcut
                vec[272] = 0;

//            if (fDimBias.state()!=BIAS::State::kRamping)
//            {
                DimClient::sendCommandNB("BIAS_CONTROL/SET_ALL_CHANNELS_VOLTAGE",
                                         vec.data(), BIAS::kNumChannels*sizeof(float));

                UdrpAvg /= 320;
                UdrpRms /= 320;
                UdrpRms -= UdrpAvg*UdrpAvg;
                UdrpRms  = UdrpRms<0 ? 0 : sqrt(UdrpRms);

                ostringstream msg;
                msg << fixed;
                msg << setprecision(2) << "dU(" << fTemp << "degC)="
                    << setprecision(3) << fTempOffsetAvg << "V+-" << fTempOffsetRms << "  Udrp="
                    << UdrpAvg << "V+-" << UdrpRms;
                msg.unsetf(ios_base::floatfield);

                if (fVoltageReduction==0)
                    msg << " Unom=" << voltageoffset << "V";
                else
                    msg << " Ured=" << fVoltageReduction << "V";

                msg << " Uov=" << Uov;
                msg << " Imed=" << data.Imed << "uA [N=" << Ndev[0] << "/" << Ndev[1] << "/" << Ndev[2] << "]";
                Info(msg);
            }
        }
        else
        {
            if (fDimBias.state()==BIAS::State::kVoltageOn)
            {
                ostringstream msg;
                msg << setprecision(4) << "Current status: dU(" << fTemp << "degC)=" << fTempOffsetAvg << "V+-" << fTempOffsetRms << ", Unom=" << voltageoffset << "V, Uov=" << (num[0]+num[1]>0?(avg[0]+avg[1])/(num[0]+num[1]):0) << " [N=" << Ndev[0] << "/" << Ndev[1] << "/" << Ndev[2] << "]";
                Info(msg);
            }
        }

        //if (GetCurrentState()>=Feedback::State::kOnStandby &&
        //    fDimBias.state()==BIAS::State::kRamping)
        //    return newstate;

        // --------------------------------- Console out --------------------------------------

        if (fIsVerbose && fDimBias.state()!=BIAS::State::kRamping)
        {
            sort(med[0].begin(), med[0].begin()+num[0]);
            sort(med[1].begin(), med[1].begin()+num[1]);

            ostringstream msg;
            msg << "   Avg0=" << setw(7) << avg[0]/num[0]    << "  |  Avg1=" << setw(7) << avg[1]/num[1];
            Debug(msg);

            msg.str("");
            msg << "   Med0=" << setw(7) << med[0][num[0]/2] << "  |  Med1=" << setw(7) << med[1][num[1]/2];
            Debug(msg);

            msg.str("");
            msg << "   Min0=" << setw(7) << min[0]           << "  |  Min1=" << setw(7) << min[1];
            Debug(msg);

            msg.str("");
            msg << "   Max0=" << setw(7) << max[0]           << "  |  Max1=" << setw(7) << max[1];
            Debug(msg);
        }

        // ---------------------------- Calibrated Currents -----------------------------------

        // FIXME:
        //  + Current overvoltage
        //  + Temp offset
        //  + User offset
        //  + Command overvoltage
        fDimCurrents.setQuality(GetCurrentState());
        fDimCurrents.setData(&data, sizeof(dim_data));
        fDimCurrents.Update(evt.GetTime());

        // FIXME: To be checked
        return GetCurrentState()==Feedback::State::kCalibrated ? Feedback::State::kCalibrated : newstate;
    }

Here is the call graph for this function:

Here is the caller graph for this function: