Index: trunk/Mars/mcore/DrsCalib.h =================================================================== --- trunk/Mars/mcore/DrsCalib.h (revision 14190) +++ trunk/Mars/mcore/DrsCalib.h (revision 14197) @@ -6,4 +6,9 @@ #include "fits.h" +#include "ofits.h" + +#ifdef __MARS__ +#include "MTime.h" +#endif class DrsCalibrate @@ -212,8 +217,14 @@ 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,14 +234,22 @@ 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,5 +274,29 @@ double pos; uint16_t cnt; + + static bool sort(const Step &s, const Step &r) { return s.avg &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,6 +304,10 @@ const int16_t offset, const uint16_t *map) { - Step rc; - + + std::vector list; + + // Fill steps into array + // Exclude broken pixels? + // Remove maximum and minimum patches (4max and 4min)? for (size_t ch=0; ch5) + { + partial_sort(list.begin(), list.begin()+list.size()/2, + list.end(), Step::sort); + + rc = AverageSteps(list, list.size()/2); + } for (size_t ch=0; chroi) + return; + + for (float *pix=vec; pix &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) {