Changeset 14197

Timestamp:

06/21/12 11:58:54 (12 years ago)

Author:

tbretz

Message:

Added code to correct for the drs steps; added WriteFitsImp and code to process a sliding average.

File:

• trunk/Mars/mcore/DrsCalib.h (modified) (8 diffs)

trunk/Mars/mcore/DrsCalib.h

@@ -6 +6 @@
 
 #include "fits.h"
+#include "ofits.h"
+
+#ifdef __MARS__
+#include "MTime.h"
+#endif
 
 class DrsCalibrate
@@ -212 +217 @@
     static double FindStep(const size_t ch0, const float *vec, int16_t roi, const int16_t pos, const uint16_t *map)
     {
-        if (pos<=1 || pos>=roi)
+        // We have about 1% of all cases which are not ahndled here,
+        // because the baseline jumps up just before the readout window
+        // and down just after it. In this cases we could determine the jump
+        // from the board time difference or throw the event away.
+        if (pos==0 || pos>=roi)
             return 0;
 
         double step = 0; // before
+        double rms  = 0; // before
+        int    cnt  = 0;
 
         // Exclude TM channel
@@ -223 +234 @@
             const size_t sw = map[hw]*roi + pos;
 
-            step += vec[sw]-vec[sw-1];
-        }
-
-        return step/8;
-    }
-
-    static void SubtractStep(const size_t ch0, const double avg, float *vec, int16_t roi, uint32_t dist, const uint16_t *map)
-    {
-        const int begin  = avg>0 ? dist :    0;
-        const int end    = avg>0 ?  roi : dist;
+            const double diff = vec[sw]-vec[sw-1];
+
+            step += diff;
+            rms  += (vec[sw]-vec[sw-1])*(vec[sw]-vec[sw-1]);
+
+            cnt++;
+        }
+
+        return cnt==0 ? 0 : step/cnt;
+    }
+
+    static void SubtractStep(const size_t ch0, const double avg, float *vec, int16_t roi, int32_t pos, const uint16_t *map)
+    {
+        if (pos==0 || pos>=roi)
+            return;
+
+        const int begin = avg>0 ? pos :   0;
+        const int end   = avg>0 ? roi : pos;
 
         const double sub = fabs(avg);
@@ -255 +274 @@
         double   pos;
         uint16_t cnt;
+
+        static bool sort(const Step &s, const Step &r) { return s.avg<r.avg; }
     };
+
+    static Step AverageSteps(const std::vector<Step> &list, uint32_t n)
+    {
+        Step rc;
+        for (auto it=list.begin(); it!=list.begin()+n; it++)
+        {
+            rc.pos += it->pos;
+            rc.avg += it->avg;
+            rc.rms += it->avg*it->avg;
+        }
+
+        rc.cnt = n;
+
+        rc.pos /= rc.cnt;
+        rc.avg /= rc.cnt;
+        rc.rms /= rc.cnt;
+
+        rc.rms = sqrt(rc.rms-rc.avg*rc.avg);
+
+        return rc;
+    }
+
 
     static Step CorrectStep(float *vec, uint16_t nch, uint16_t roi,
@@ -261 +304 @@
                             const int16_t offset, const uint16_t *map)
     {
-        Step rc;
-
+
+        std::vector<Step> list;
+
+        // Fill steps into array
+        // Exclude broken pixels?
+        // Remove maximum and minimum patches (4max and 4min)?
         for (size_t ch=0; ch<nch; ch += 9)
         {
@@ -273 +320 @@
                 continue;
 
-            rc.pos += dist;
-            rc.avg += step;
-            rc.rms += step*step;
-            rc.cnt++;
-        }
-
-        if (rc.cnt==0 || rc.avg==0)
+            Step rc;
+            rc.pos = dist;
+            rc.avg = step;
+            list.push_back(rc);
+        }
+
+        if (list.size()==0)
             return Step();
 
-        rc.pos /= rc.cnt;
-        rc.avg /= rc.cnt;
-        rc.rms /= rc.cnt;
-
-        rc.rms = sqrt(rc.rms-rc.avg*rc.avg);
+        Step rc = AverageSteps(list, list.size());;
+
+        if (rc.avg==0)
+            return Step();
+
+        if (rc.rms>5)
+        {
+            partial_sort(list.begin(), list.begin()+list.size()/2,
+                         list.end(), Step::sort);
+
+            rc = AverageSteps(list, list.size()/2);
+        }
 
         for (size_t ch=0; ch<nch; ch += 9)
@@ -367 +421 @@
                     i += 2;
                 }
+            }
+        }
+    }
+
+    static void SlidingAverage(float *const vec, const uint32_t roi, const uint16_t w)
+    {
+        if (w==0)
+            return;
+
+        if (w>roi)
+            return;
+
+        for (float *pix=vec; pix<vec+1440*roi; pix += roi)
+        {
+            for (float *ptr=pix; ptr<pix+roi-w; ptr++)
+            {
+                for (float *p=ptr+1; p<ptr+w; p++)
+                    *ptr += *p;
+                *ptr /= w;
             }
         }
@@ -884 +957 @@
     }
 
+    std::string WriteFitsImp(const std::string &filename, const std::vector<float> &vec) const
+    {
+        const size_t n = 1440*1024*4 + 1440*fRoi*2 + fNumTm*fRoi*2 + 3;
+
+        std::ofits file(filename.c_str());
+        if (!file)
+        {
+            std::ostringstream msg;
+            msg << "Could not open file '" << filename << "': " << strerror(errno);
+            return msg.str();
+        }
+
+        file.AddColumnInt("RunNumberBaseline");
+        file.AddColumnInt("RunNumberGain");
+        file.AddColumnInt("RunNumberTriggerOffset");
+
+        file.AddColumnFloat(1024*1440,   "BaselineMean",        "mV");
+        file.AddColumnFloat(1024*1440,   "BaselineRms",         "mV");
+        file.AddColumnFloat(1024*1440,   "GainMean",            "mV");
+        file.AddColumnFloat(1024*1440,   "GainRms",             "mV");
+        file.AddColumnFloat(fRoi*1440,   "TriggerOffsetMean",   "mV");
+        file.AddColumnFloat(fRoi*1440,   "TriggerOffsetRms",    "mV");
+        file.AddColumnFloat(fRoi*fNumTm, "TriggerOffsetTMMean", "mV");
+        file.AddColumnFloat(fRoi*fNumTm, "TriggerOffsetTMRms",  "mV");
+
+#ifdef __MARS__
+        const MTime now;
+        file.SetStr(  "TELESCOP", "FACT",               "Telescope that acquired this data");
+        file.SetStr(  "PACKAGE",  "MARS",               "Package name");
+        file.SetStr(  "VERSION",  "1.0",                "Package description");
+        //file.SetStr(  "CREATOR",  "root",               "Program that wrote this file");
+        file.SetFloat("EXTREL",   1.0,                  "Release Number");
+        file.SetStr(  "COMPILED", __DATE__" "__TIME__,  "Compile time");
+        //file.SetStr(  "REVISION", REVISION,             "SVN revision");
+        file.SetStr(  "ORIGIN",   "FACT",               "Institution that wrote the file");
+        file.SetStr(  "DATE",     now.GetStringFmt("%Y-%m-%dT%H:%M:%S").Data(), "File creation date");
+        file.SetInt(  "NIGHT",    now.GetNightAsInt(),  "Night as int");
+        file.SetStr(  "TIMESYS",  "UTC",                "Time system");
+        file.SetStr(  "TIMEUNIT", "d",                  "Time given in days w.r.t. to MJDREF");
+        file.SetInt(  "MJDREF",   40587,                "MJD to UNIX time (seconds since 1970/1/1)");
+#else
+        DataWriteFits2::WriteDefaultKeys(file);
+#endif
+        file.SetInt("STEP", fStep, "");
+
+        file.SetInt("ADCRANGE", 2000,    "Dynamic range of the ADC in mV");
+        file.SetInt("DACRANGE", 2500,    "Dynamic range of the DAC in mV");
+        file.SetInt("ADC",      12,      "Resolution of ADC in bits");
+        file.SetInt("DAC",      16,      "Resolution of DAC in bits");
+        file.SetInt("NPIX",     1440,    "Number of channels in the camera");
+        file.SetInt("NTM",      fNumTm,  "Number of time marker channels");
+        file.SetInt("NROI",     fRoi,    "Region of interest");
+
+        file.SetInt("NBOFFSET", fNumOffset,       "Num of entries for offset calibration");
+        file.SetInt("NBGAIN",   fNumGain/1953125, "Num of entries for gain calibration");
+        file.SetInt("NBTRGOFF", fNumTrgOff,       "Num of entries for trigger offset calibration");
+
+        // file.WriteKeyNT("DAC_A",    fData.fDAC[0],    "Level of DAC 0 in DAC counts")   ||
+        // file.WriteKeyNT("DAC_B",    fData.fDAC[1],    "Leval of DAC 1-3 in DAC counts") ||
+        // file.WriteKeyNT("DAC_C",    fData.fDAC[4],    "Leval of DAC 4-7 in DAC counts") ||
+
+        file.WriteTableHeader("DrsCalibration");
+
+        if (!file.WriteRow(vec.data()+1, n*sizeof(float)))
+        {
+            std::ostringstream msg;
+            msg << "Writing data to " << filename << " failed.";
+            return msg.str();
+        }
+
+        return std::string();
+    }
+
+
     std::string ReadFitsImp(const std::string &str)
     {