Python으로 Ortec Chn 파일 읽기

 

Ortec 검출기로 측정한 스펙트럼 CHN 파일을 Python으로 읽어내기

 

이 글을 쓰려고 보니, 정말 특정 분야에만 이용될테고 과연 이 코드를 사용할 사람이 전국에 몇명이나 될지 모를만한 아무 쓸모없는 글을 올리는 듯한 생각이 듭니다. 그래도 언젠가는 이게 도움이 되는 분이 있겠죠.

 

방사선 중 감마선을 측정하여 핵종 분석하는 분들은 HPGe 검출기들을 주로 쓰실테고, 측정 목적 및 용도에 따라 NaI와 같은 신틸레이터 검출기를 쓰기도 합니다. 쓰고보니, 비단 감마선 뿐만 아니라 알파선이나 베타선도 목적에 따라서는 채널(에너지)에 따른 카운트가 얼마나 쌓였는지 히스토그램의 형태로 스펙스럼을 봐야할 일이 생깁니다.

 

전문가 분들이야 알아서 필요한 장비들을 사고, DAQ와 관련한 코드를 직접 만드시겠지만, 감마선 핵종 분석과 같은 분야에서는 굳이 그러지 않아도 잘 만들어진 프로그램들이 있습니다. 어차피 흔히 이용되는 검출기는 ORTEC사나 CANBERRA사의 검출기를 쓸테고, 제공되는 프로그램을 쓰면 되니까요. 요즘에는 BSI사의 제품도 거론되는 것 같은데, 실제로 사용하는 곳은 본 적이 없습니다.

 

Ortec사의 검출기를 이용하시는 분들은 아래 그림과 같은 프로그램을 많이 봤을 겁니다. ORTEC사의 GammaVision 프로그램의 모습입니다.

GammaVision 프로그램(출처: https://www.ortec-online.com)

 

Canberra사의 검출기는 제가 다뤄보지는 않아서 모르겠습니다만, Ortec사의 HPGe검출기를 사용하면서 Maestro프로그램으로 데이터를 받아서 GammaVision으로 분석하는 것을 지긋지긋하게 많이 구경해봤습니다. (정작 저는 GammaVision을 가끔 다루기만 했습니다..)

 

GammaVison이 참 좋은 프로그램이고 기능도 좋습니다만, 늘 정해진 핵종의 peak만 볼 것이 아니고, 정말 다양한 상황에 놓이게 되다보면 좀 더 자유도가 높은 무언가가 필요하게 됩니다. 뭐 거창한 것은 아니고, 그냥 직접 코드짜서 데이터를 읽고, 직접 데이터로 지지고 볶는 것이죠.

 

랩미팅에서 늘 보던 GammaVision에서 뽑아준 스펙트럼 그래프도 질리고, 뭔가 프로그램에 아쉬운 점이 있고 (2,3개 peak이 겹쳤을때, 피팅하는 방법이 요즘 버전은 어떤지 모르겠지만 제가 사용하던 버전에서는...참..), 어차피 그래프는 CERN ROOT로 그리던 시절이라 데이터 분석을 pyROOT로 할겸 CHN파일을 읽는 코드를 파이썬으로 만들었었습니다.

 

Github에 올려서 다운로드 받게하면 더 좋겠지만, 제가 공개된 git을 쓰지 않습니다. 제 개인 코드는 집안에서만 개인서버의 git을 쓰기 때문에 그냥 무식하게 코드를 올립니다.

 

Ortec Maestro로 생성된 CHN 파일 읽기

# Only for new version of CHN

import struct
import datetime

p_ver = "v0.001"

class OrtecChn:
  def __init__(self):
    # Begining of the CHN file
    self.INDEX_CHN = 0
    self.DET_N = 0
    self.SEG_N = 0
    self.ASCII_START_SS = ""
    self.REALTIME = 0.
    self.LIVETIME = 0.
    self.ASCII_START_DATE = ""
    self.ASCII_START_HHMM = ""
    self.CHN_OFFSET = 0
    self.CHN_N = 0
    self.DATA = []
 
    # New Version
    self.INDEX_NEW = 0
    self.RESERV2 = ""
    self.E_CAL_0_INTERCEPT = 0.
    self.E_CAL_SLOPE = 1.
    self.E_CAL_QUAD_TERM = 0.
    self.P_SHAPE_CAL_0_INTERCEPT = 1.
    self.P_SHAPE_CAL_SLOPE = 0.
    self.P_SHAPE_CAL_QUAD_TETM = 0.
    self.RESERV28 = ""
    self.LEN_DET_DSC = 0
    self.DET_DSC = ""
    self.LEN_SAM_DSC = 0
    self.SAM_DSC = ""
    self.RESERV384 = ""

    # ----
    self.fChn_filename = ""
    self.fChn_input = ""
    self.fChn_outname = ""
    self.fRealtime = 0.
    self.fLivetime = 0.
    self.fDead = 0.
    self.ascii_date = ""
    self.ascii_time = ""

  def OpenChn (self, filename):
    self.fChn_filename = filename
    self.fChn_input = open(self.fChn_filename, "rb")
    print ("Open "+self.fChn_filename)
    self.ReadChn()

  def GetData(self):
    return self.DATA

  def ReadChn (self):
    #----- 1st part
    self.INDEX_CHN = int(struct.unpack('h', self.fChn_input.read(2))[0])
    if not self.INDEX_CHN == -1:
      print ("This is not a Ortec Chn file.")
      exit()

    self.DET_N = int(struct.unpack('h', self.fChn_input.read(2))[0])
    self.SEG_N = int(struct.unpack('h', self.fChn_input.read(2))[0])

    for i in range (2):
      tmp_ss = str(struct.unpack('s', self.fChn_input.read(1))[0])
      self.ASCII_START_SS = self.ASCII_START_SS + tmp_ss
    
    self.REALTIME = int(struct.unpack('i', self.fChn_input.read(4))[0])
    self.fRealtime = self.ConvertTime(self.REALTIME)
    self.LIVETIME = int(struct.unpack('i', self.fChn_input.read(4))[0])
    self.fLivetime = self.ConvertTime(self.LIVETIME)
    self.fDead = (self.fRealtime - self.fLivetime)/self.fRealtime*100
    # print (self.fDead)

    for i in range (8):
      tmp_date = str(struct.unpack('s',self.fChn_input.read(1))[0]) 
      self.ASCII_START_DATE = self.ASCII_START_DATE + tmp_date
    for i in range (4):
      tmp_hhmm = str(struct.unpack('s',self.fChn_input.read(1))[0]) 
      self.ASCII_START_HHMM = self.ASCII_START_HHMM + tmp_hhmm

    self.CHN_OFFSET = int(struct.unpack('h', self.fChn_input.read(2))[0])
    self.CHN_N = int(struct.unpack('h', self.fChn_input.read(2))[0])

    for i in range (self.CHN_N):
      tmp = int(struct.unpack('i', self.fChn_input.read(4))[0])
      self.DATA.append(tmp)

    #------ 2nd part
    self.INDEX_NEW = int(struct.unpack('h', self.fChn_input.read(2))[0])
    if not self.INDEX_NEW == -102:
      print ("This CHN file is not a new version")
      exit()
    self.RESERV2 = self.fChn_input.read(2)
    self.E_CAL_0_INTERCEPT = float(struct.unpack('f', self.fChn_input.read(4))[0])
    self.E_CAL_SLOPE = float(struct.unpack('f', self.fChn_input.read(4))[0])
    self.E_CAL_QUAD_TERM = float(struct.unpack('f', self.fChn_input.read(4))[0])
    self.P_SHAPE_CAL_0_INTERCEPT = float(struct.unpack('f', self.fChn_input.read(4))[0])
    self.P_SHAPE_CAL_SLOPE = float(struct.unpack('f', self.fChn_input.read(4))[0])
    self.P_SHAPE_CAL_QUAD_TETM = float(struct.unpack('f', self.fChn_input.read(4))[0])

    self.RESERV28 = self.fChn_input.read(228)
    self.LEN_DET_DSC =  int(struct.unpack('b', self.fChn_input.read(1))[0])
    for i in range(63):
      tmp = str(struct.unpack('s', self.fChn_input.read(1))[0])
      self.DET_DSC = self.DET_DSC + tmp
    self.LEN_SAM_DSC = int (struct.unpack('b', self.fChn_input.read(1))[0])
    for i in range(63):
      tmp = str(struct.unpack('s', self.fChn_input.read(1))[0])
      self.SAM_DSC = self.SAM_DSC + tmp
      
    self.RESERV384 = self.fChn_input.read(128)

  def ConvertTime(self, ftime):
    return ftime*0.020
  
  def ConvertDate(self):
    self.ascii_date = self.ASCII_START_DATE[0:2]+"/"+self.ASCII_START_DATE[2:5]+"/"
    if self.ASCII_START_DATE[-1] == "1":
      self.ascii_date = self.ascii_date+"20"+self.ASCII_START_DATE[5:7]

  def WriteTXT(self):
    self.fChn_outname = self.fChn_filename.split(".Chn")[0] + ".txt"
    fout = open(self.fChn_outname, "w")
    fout.write("$SPEC_ID: \n")
    fout.write(self.SAM_DSC[:self.LEN_SAM_DSC]+"\n")
    fout.write("$SPEC_REM: \n")
    fout.write("DETECTOR NUMBER# "+ str(self.DET_N)+"\n")
    fout.write("DETDESC# "+self.DET_DSC[0:self.LEN_DET_DSC]+"\n")
    fout.write("AP# " + p_ver+"\n")
    fout.write("$DATE_MEASUREMENT:\n")
    self.ConvertDate()
    fout.write(self.ascii_date+"\t"+self.ASCII_START_HHMM[0:2]+":"+self.ASCII_START_HHMM[2:4]+":"+self.ASCII_START_SS+"\n")
    fout.write("$MEASUREMENT_TIME:\n")
    fout.write(str(self.fLivetime)+"\t"+str(self.fRealtime)+"\n")
    fout.write("$DATA: \n")
    fout.write(str(self.CHN_OFFSET)+"\t"+str(self.CHN_N)+"\n")
    for i in range(self.CHN_N):
      fout.write(str(self.DATA[i])+"\n")

  def CloseChn(self):
    self.fChn_input.close()
    print ("Close "+self.fChn_filename)

  def GetRealtime(self):
    return self.fRealtime 

  def GetLivetime(self):
    return self.fLivetime

  def GetMonth(self,strMonth):
    monList = ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"]
    for i in range(len(monList)):
      if strMonth == monList[i]:
        return i+1

  def GetStartTime(self):
    self.ConvertDate()
    fyear = int(self.ascii_date.split("/")[2])
    fmonth = self.GetMonth(self.ascii_date.split("/")[1])
    fday = int(self.ascii_date.split("/")[0])
    fhh = int(self.ASCII_START_HHMM[0:2])
    fmm = int(self.ASCII_START_HHMM[2:])
    fss = int(self.ASCII_START_SS)
    startTime = datetime.datetime(fyear,fmonth,fday,fhh,fmm,fss)
    return startTime

if __name__ == "__main__":
  fChn = OrtecChn()
  fChn.OpenChn("test_21600s.Chn")
  print(fChn.GetStartTime())
  print(fChn.GetLivetime())
  print(fChn.GetRealtime())
  fChn.WriteTXT()
  fChn.CloseChn()

 

수년이 지나 지금와서 코드를 보니 꽤 무식하게 코드를 짠 부분이 곳곳에 보이네요. Get함수가 몇개 없는걸 보아하니, 바로 접근해서 썼던 것 같습니다. 기억이..

NAS에 옛 자료들을 정리하다 눈에 띄었던 파일으로, 2015년도에 작성된거라 요즘에도 잘 되련지 모르겠습니다. 아마 파이썬3을 쓰긴 했었던 것 같고, Ortec의 Chn파일 형식이 그대로인지 모르겠군요. 5년사이에 바꾸진 않았겠죠.