週末はいつも晴れ

社会人3年目の日記です。プログラミングとか旅行とかラーメン。

エネルギー較正時にカウントをバンチングするクラス[C++]

タイトル通りのことをやりたいのですが、絶対意味が伝わらないorz
日本語って難しいですね。


具体的には次のような例です:
f:id:ikarino99:20140705124058p:plain
一般的に放射線などに起因する電圧がデータ収集系に取り込まれる際には、電圧に比例した離散化したデジタルチャンネルに入力されます。その後、実験前(後)に標準線源などを用いてチャンネルに対応するエネルギーの較正を行います。
つまり、こういうことです:
f:id:ikarino99:20140705125203p:plain
一般的にエネルギー較正時のx軸の値は、区切れのよい数値ではありません。

その後の解析や、グラフ化の際にいびつな形になってしまうのを防ぐために、x軸を変換する必要があります。



区切れのよい数値に変換する方法はいくつかあります。

一つの方法はスプライン補間を行う方法です。
スプライン補間とは、得られたデータ点をそれぞれの区間で異なる3次関数で結び、点と点の間を補間する方法です。
gnuplotの with line コマンドではスプライン補間が行われます。
スプライン補間した後、必要な区切りの良いxにおけるyの値を読み取れば良いのです。

非常にわかりやすい方法ですが、スプライン補間には得られた統計量を反映することができないという弱点があります。
特殊なものをのぞき、放射線測定では統計量を起源とする誤差がほとんどです。
誤差を小さくするには統計量を無駄にすることはできないのです。
また、スプライン補間では細かな変化に敏感になるため、変化の激しいスペクトル間隔を広げる変換の場合、全く正当性が保証されなくなります。



区切れのよい数値に変換しつつ、得られた統計量を無駄にしないバンチングを行うには、簡単に言うと足し算を行う必要があります。

100keV刻みにデータを変換したいとすると、150keV〜250keVに含まれる元のカウントを200keVのカウントとして和をとるということです。

こう説明すると簡単に感じますが、複数の領域に重なってしまった場合に分配する際に割と面倒な計算をする必要があります。

ROOTにHistgramクラスのメンバ関数として準備されていればよかったのですが、見つからなかったので仕方なく自分で書くことにしました。

// Copyright 2014 ikarino
/*
 * ver2.cpp: エネルギー較正時にカウントをバンチングするC++プログラム
 *           C++11を使っているので、コンパイル時は-std=c++11オプションを。
 */
#include <cstdio>
#include <string>
#include <vector>
#include <fstream>
// #define DEBUG

/* class: Data
 * チャンネル(エネルギー)とカウントを1つずつ持っているクラス。
 */
class Data {
 public:
  Data(double channel, double count)
      : _channel(channel),
        _count(count)
  {}
  void print() {
    printf("%f\t%f\n", _channel, _count);
  }
  double getChannel() {return _channel;}
  double getCount() {return _count;}
 private:
  double _channel;
  double _count;
};

/* class: Hist
 * Dataを複数格納するクラス
 */
class Hist {
 public:
  Hist() {
    dataSet = new std::vector<Data>;
  }
  explicit Hist(std::string filename) {
    dataSet = new std::vector<Data>;
    this->ReadFile(filename);
  };
  ~Hist() {
    delete dataSet;
  }
  bool ResetBin(double min, double max, int bin_num);
  void ReadFile(std::string filename);
  void SaveFile(std::string filename);
 private:
  std::vector<Data> *dataSet;
};

void Hist::SaveFile(std::string filename) {
  std::ofstream ofs(filename, std::ios::out);
  for (auto data = dataSet->begin(); data != dataSet->end(); data++) {
    ofs << data->getChannel() << "\t" << data->getCount() << "\n";
  }
}

void Hist::ReadFile(std::string filename) {
  std::ifstream ifs(filename);
  double x, y;
  while (ifs >> x >> y) {
    Data s(x, y);
    dataSet->push_back(s);
  }
}

/*
 * 主要部分
 * 新しいエネルギーの最小値、最大値、データ点数を受け取り、新しいスペクトルをdataSetに入れる。
 */
bool Hist::ResetBin(double min, double max, int bin_num) {
  if (bin_num < 2) {
    printf("bin number must be > 2\n");
    return false;
  }
  if (bin_num > 65536) {
    return false;
  }

  std::vector<double> breakline(bin_num+1);
  double bin_width = (max-min)/(bin_num-1);
  for (int i = 0; i <= bin_num; i++) {
    breakline.at(i) = min - bin_width/2 + bin_width*i;
  }
#ifdef DEBUG
  printf("bin_width = %f\n", bin_width);
#endif
  std::vector<double> bunchedCount(bin_num, 0);

  for (auto data = dataSet->begin(); data != dataSet->end(); data++) {
    /*
     * x軸の下限と上限の処理
     * 線形的に決定しているので、元のx軸が線形的でない場合は注意。
     *    |   |   |   |   |   |
     *  : | : | : | : | : | : | :
     *  a 1 b 2 b 3 b 4 b 5 b 6 c
     *  
     *  b: 中間値を採用
     *  a: 1と2の差の半分を1から引いた。cも同様。
     * 
     */
    double bin_min = (data->getChannel() + (data-1)->getChannel())/2;
    double bin_max = (data->getChannel() + (data+1)->getChannel())/2;
    if (data == dataSet->begin()) {
      bin_min = (3 * data->getChannel() - (data+1)->getChannel())/2;
    } else if (data == dataSet->end()-1) {
      bin_max = (3 * data->getChannel() - (data-1)->getChannel())/2;
    }
    /*
     * bin_min と bin_max がbreaklineのどこにいるのか、判定を行い、placeに代入。
     *               |      |     |     |     |     |  ... |
     *            :  |   :  |  :  |  :  |  :  |  :  |  ... |  :
     * breakline  0      1     2     3     4     5         bin_num
     *          ^   ^      ^          ^                          ^
     * place   -1   0      1          3                         -2
     */
    int place_min = -1, place_max = -1;
    for (int i = 0; i <= bin_num; i++) {
      if (bin_min > breakline.at(i)) {
        place_min = i;
        if (place_min == bin_num) place_min = -2;
      }
    }
    for (int i = 0; i <=bin_num; i++) {
      if (bin_max > breakline.at(i)) {
        place_max = i;
        if (place_max == bin_num) place_max = -2;
      }
    }
#ifdef DEBUG
    printf("min %d\tmax %d\n", place_min, place_max);
#endif
    /*
     * placeにあわせてbunchedCountにカウントを加えていく。
     */
    double unitCount = data->getCount()/(bin_max - bin_min);
    if (place_min >= 0 && place_max >= 0) {
      if (place_min != place_max) {
        bunchedCount.at(place_min)
            += (breakline.at(place_min+1)-bin_min)*unitCount;
        for (int i = place_min+1; i < place_max; i++) {
          bunchedCount.at(i) += bin_width * unitCount;
        }
        bunchedCount.at(place_max)
            += (bin_max - breakline.at(place_max))*unitCount;
      } else {
        bunchedCount.at(place_min) += data->getCount();
      }
    } else if (place_min == -2 || place_max == -1) {
      // Range over なので何もしない。
    } else if (place_min != -1 && place_max == -2) {
      /*
       * こういう感じのとき。
       *               |     |  ... |
       *            :  |  :  |  ... |  :
       * breakline  4     5         bin_num
       *                             ^    ^     
       * place                      min  max
       */
      bunchedCount.at(place_min)
          += (breakline.at(place_min+1)-bin_min)*unitCount;
      for (int i = place_min+1; i < bin_num; i++) {
        bunchedCount.at(i) += bin_width * unitCount;
      }
    } else if (place_max != -2 && place_min == -1) {
      // 上の逆ver
      for (int i = 0; i < place_max; i++) {
        bunchedCount.at(i) += bin_width * unitCount;
      }
      bunchedCount.at(place_max)
          += (bin_max - breakline.at(place_max))*unitCount;
    } else {
      printf("Exception occured!\n");
      exit(1);
    }
#ifdef DEBUG
    for (int i = 0; i < bin_num; i++) {
      printf("%d\t%f\n", i, bunchedCount.at(i));
    }
#endif
  }
  dataSet->clear();
  for (int i = 0; i < bin_num; i++) {
    Data d(min + bin_width*i, bunchedCount.at(i));
    dataSet->push_back(d);
  }
  return true;
}



int main(int argc, char *argv[]) {
  Hist *dataSet = new Hist();
  dataSet->ReadFile("data");
  // dataSet->ResetBin(0, 55, 56);  // OK
  dataSet->ResetBin(0, 60, 3);  // OK
  dataSet->ResetBin(10, 30, 8);
  dataSet->SaveFile("out");
  return 0;
}