Index: trunk/FACT++/src/ratecontrol.cc
===================================================================
--- trunk/FACT++/src/ratecontrol.cc	(revision 15064)
+++ trunk/FACT++/src/ratecontrol.cc	(revision 15069)
@@ -8,4 +8,5 @@
 #include "Configuration.h"
 #include "Console.h"
+#include "PixelMap.h"
 
 #include "tools.h"
@@ -45,4 +46,6 @@
     map<string, config> fRunTypes;
 
+    PixelMap fMap;
+
     bool fTriggerOn;
 
@@ -65,5 +68,7 @@
     uint16_t fRequiredEvents;
 
-    deque<pair<Time,float>> fCurrents;
+    deque<pair<Time,float>> fCurrentsMed;
+    deque<pair<Time,float>> fCurrentsDev;
+    deque<pair<Time,vector<float>>> fCurrentsVec;
 
     bool fVerbose;
@@ -404,13 +409,19 @@
         const Time &time = evt.GetTime();
         const float med  = evt.Get<float>(416*4+4+4);
+        const float dev  = evt.Get<float>(416*4+4+4+4);
+        const float *cur = evt.Ptr<float>();
 
         // Keep all median currents of the past 10 seconds
-        fCurrents.push_back(make_pair(time, med));
-        while (!fCurrents.empty())
-        {
-            if (time-fCurrents.front().first<boost::posix_time::seconds(fAverageTime))
+        fCurrentsMed.push_back(make_pair(time, med));
+        fCurrentsDev.push_back(make_pair(time, dev));
+        fCurrentsVec.push_back(make_pair(time, vector<float>(cur, cur+320)));
+        while (!fCurrentsMed.empty())
+        {
+            if (time-fCurrentsMed.front().first<boost::posix_time::seconds(fAverageTime))
                 break;
 
-            fCurrents.pop_front();
+            fCurrentsMed.pop_front();
+            fCurrentsDev.pop_front();
+            fCurrentsVec.pop_front();
         }
 
@@ -423,5 +434,5 @@
 
         // We want at least 8 values for averaging
-        if (fCurrents.size()<fRequiredEvents)
+        if (fCurrentsMed.size()<fRequiredEvents)
             return GetCurrentState();
 
@@ -429,16 +440,74 @@
         double avg = 0;
         double rms = 0;
-        for (auto it=fCurrents.begin(); it!=fCurrents.end(); it++)
+        for (auto it=fCurrentsMed.begin(); it!=fCurrentsMed.end(); it++)
         {
             avg += it->second;
             rms += it->second*it->second;
         }
-        avg /= fCurrents.size();
-        rms /= fCurrents.size();
-
+        avg /= fCurrentsMed.size();
+        rms /= fCurrentsMed.size();
         rms = sqrt(rms-avg*avg);
 
-        fThresholdMin = max(uint16_t(36.0833*pow(avg, 0.638393)+184.037), fThresholdReference);
+        double avg_dev = 0;
+        for (auto it=fCurrentsDev.begin(); it!=fCurrentsDev.end(); it++)
+            avg_dev += it->second;
+        avg_dev /= fCurrentsMed.size();
+
+        // One could recalculate the median of all pixels incluing the
+        // correction for the three crazy pixels, but that is three out
+        // of 320. The effect on the median should be negligible anyhow.
+        vector<double> vec(160);
+        for (auto it=fCurrentsVec.begin(); it!=fCurrentsVec.end(); it++)
+            for (int i=0; i<320; i++)
+            {
+                const PixelMapEntry &hv = fMap.hv(i);
+                if (!hv)
+                    continue;
+
+                // The current is proportional to the rate. To calculate
+                // a measure for the rate, the average current per pixel
+                // is caluclated for the trigger patch.
+                int weight = hv.group() ? 5 : 4;
+
+                // Use only the current in the pixels with the correct
+                // resistor as a reference, ignore the crazy ones.
+                // Effects of these should be corrected by the
+                // rate control later, not the initial setup.
+                if (i==66)           
+                    weight = 4./(3+10);
+                if (i==191 || i==193)
+                    weight = 5./(4+10);
+
+                vec[hv.hw()/9] += it->second[i] * weight;
+            }
+
+        //fThresholdMin = max(uint16_t(36.0833*pow(avg, 0.638393)+184.037), fThresholdReference);
+        fThresholdMin = max(uint16_t(36.0833*pow(avg, 0.638393)+210), fThresholdReference);
         fThresholds.assign(160, fThresholdMin);
+
+        const int32_t val[2] = { -1, fThresholdMin };
+        Dim::SendCommand("FTM_CONTROL/SET_THRESHOLD", val);
+
+        double avg2 = 0;
+        for (int i=0; i<160; i++)
+        {
+            vec[i] /= fCurrentsVec.size()*9;
+            avg2 += vec[i];
+
+            if (vec[i]-avg>6*avg_dev)
+            {
+                fThresholds[i] = max(uint16_t(36.0833*pow(vec[i], 0.638393)+185), fThresholdReference);
+
+                cout << "i=" << i << " " << avg << " " << avg_dev << " " << vec[i] << " " << fThresholds[i] << endl;
+
+                const int32_t dat[2] = { i, fThresholds[i] };
+                Dim::SendCommand("FTM_CONTROL/SET_THRESHOLD", dat);
+
+                fBlock[i/4] = true;
+            }
+        }
+
+        avg2 /= 160;
+        cout << "AVG2="<<avg2 << endl;
 
         const RateControl::DimThreshold data = { fThresholdMin, fCalibrationTimeStart.Mjd(), Time().Mjd() };
@@ -448,5 +517,5 @@
         ostringstream out;
         out << setprecision(3);
-        out << "Measured average current " << avg << "uA +- " << rms << "uA [N=" << fCurrents.size() << "]... mininum threshold set to " << fThresholdMin;
+        out << "Measured average current " << avg << "uA +- " << rms << "uA [N=" << fCurrentsMed.size() << "]... mininum threshold set to " << fThresholdMin;
         Info(out);
 
@@ -470,4 +539,5 @@
         fTriggerRate  = -1;
         fCounter      = 0;
+        fBlock.assign(160, false);
 
         fCalibrateByCurrent = false;
@@ -495,7 +565,8 @@
         fCalibrateByCurrent = true;
         fCalibrationTimeStart = Time();
+        fBlock.assign(160, false);
 
         ostringstream out;
-        out << "Rate calibration by current " << fThresholdReference << " with a target rate of " << fTargetRate << " Hz";
+        out << "Rate calibration by current " << fThresholdReference << ".";
         Info(out);
 
@@ -507,9 +578,15 @@
         const string name = evt.GetText();
 
-        const auto it = fRunTypes.find(name);
+        auto it = fRunTypes.find(name);
         if (it==fRunTypes.end())
         {
-            Error("CalibrateRun - Run-type '"+name+"' not found.");
-            return GetCurrentState();
+            Info("CalibrateRun - Run-type '"+name+"' not found... trying 'default'.");
+
+            it = fRunTypes.find("default");
+            if (it==fRunTypes.end())
+            {
+                Error("CalibrateRun - Run-type 'default' not found.");
+                return GetCurrentState();
+            }
         }
 
@@ -698,8 +775,11 @@
     }
 
-    bool CheckConfig(Configuration &conf, const string &name, const string &sub)
+    bool GetConfig(Configuration &conf, const string &name, const string &sub, uint16_t &rc)
     {
         if (conf.HasDef(name, sub))
+        {
+            rc = conf.GetDef<uint16_t>(name, sub);
             return true;
+        }
 
         Error("Neither "+name+"default nor "+name+sub+" found.");
@@ -710,4 +790,10 @@
     {
         fVerbose = !conf.Get<bool>("quiet");
+
+        if (!fMap.Read(conf.Get<string>("pixel-map-file")))
+        {
+            Error("Reading mapping table from "+conf.Get<string>("pixel-map-file")+" failed.");
+            return 1;
+        }
 
         fThresholdReference = 300;
@@ -735,37 +821,11 @@
             }
 
-            if (!CheckConfig(conf, "calibration-type.", *it))
+            config &c = fRunTypes[*it];
+            if (!GetConfig(conf, "calibration-type.", *it, c.fCalibrationType) ||
+                !GetConfig(conf, "target-rate.",      *it, c.fTargetRate)      ||
+                !GetConfig(conf, "min-threshold.",    *it, c.fMinThreshold)    ||
+                !GetConfig(conf, "average-time.",     *it, c.fAverageTime)     ||
+                !GetConfig(conf, "required-events.",  *it, c.fRequiredEvents))
                 return 2;
-
-            config c;
-            c.fCalibrationType = conf.GetDef<uint16_t>("calibration-type.",  *it);
-
-            switch (c.fCalibrationType)
-            {
-            case 0:
-                break;
-
-                // Calibrate by rate
-            case 1:
-                if (!CheckConfig(conf, "target-rate.",   *it) ||
-                    !CheckConfig(conf, "min-threshold.", *it))
-                    return 3;
-                c.fTargetRate   = conf.GetDef<uint16_t>("target-rate.",    *it);
-                c.fMinThreshold = conf.GetDef<uint16_t>("min-threshold.", *it);
-                break;
-
-                // Calibrate by current
-            case 2:
-                if (!CheckConfig(conf, "min-threshold.",   *it) ||
-                    !CheckConfig(conf, "average-time.",    *it) ||
-                    !CheckConfig(conf, "required-events.", *it))
-                    return 4;
-                c.fMinThreshold    = conf.GetDef<uint16_t>("min-threshold.",   *it);
-                c.fAverageTime     = conf.GetDef<uint16_t>("average-time.",    *it);
-                c.fRequiredEvents  = conf.GetDef<uint16_t>("required-events.", *it);
-                break;
-            }
-
-            fRunTypes[*it] = c;
         }
 
@@ -788,5 +848,6 @@
     po::options_description control("Rate control options");
     control.add_options()
-        ("quiet,q",       po_bool(),  "Disable printing more informations during rate control.")
+        ("quiet,q", po_bool(),  "Disable printing more informations during rate control.")
+        ("pixel-map-file", var<string>("FACTmapV5a.txt"), "Pixel mapping file. Used here to get the default reference voltage.")
        //("max-wait",   var<uint16_t>(150), "The maximum number of seconds to wait to get the anticipated resolution for a point.")
        // ("resolution", var<double>(0.05) , "The minimum resolution required for a single data point.")
@@ -797,10 +858,10 @@
     po::options_description runtype("Run type configuration");
     runtype.add_options()
-        ("run-type",           vars<string>(),   "Name of run-types (replace the * in the following configuration by the case-sensitive names defined here)")
-        ("calibration-type.*", vars<uint16_t>(), "Calibration type (0: none, 1: by rate, 2: by current)")
-        ("target-rate.*",      vars<uint16_t>(), "Target rate for calibration by rate")
-        ("min-threshold.*",    vars<uint16_t>(), "Minimum threshold which can be applied in a calibration")
-        ("average-time.*",     vars<uint16_t>(), "Time in seconds to average the currents for a calibration by current.")
-        ("required-events.*",  vars<uint16_t>(), "Number of required current events to start a calibration by current.");
+        ("run-type",           vars<string>(),  "Name of run-types (replace the * in the following configuration by the case-sensitive names defined here)")
+        ("calibration-type.*", var<uint16_t>(), "Calibration type (0: none, 1: by rate, 2: by current)")
+        ("target-rate.*",      var<uint16_t>(), "Target rate for calibration by rate")
+        ("min-threshold.*",    var<uint16_t>(), "Minimum threshold which can be applied in a calibration")
+        ("average-time.*",     var<uint16_t>(), "Time in seconds to average the currents for a calibration by current.")
+        ("required-events.*",  var<uint16_t>(), "Number of required current events to start a calibration by current.");
     ;