tools.py

import math
import numpy as np
import json
import re
from settings import *

def PMT_setup(pmt_file_with_index):
  '''
  Reads in PMT positional information according to PMT location file.
  This file is internal to KamLAND collaboration thus cannot be provided here,
  it should follow the format of:

  PMTID  PMT_X[cm]  PMT_Y[cm]  PMT_Z[cm]

  Each line represents a PMT and each field is separated by blank space
  '''
  PMT_POSITION = {}
  for pmt in np.loadtxt(pmt_file_with_index):
    current_pmt_pos = pmt[1:] / 100.0
    PMT_POSITION[int(pmt[0])] = current_pmt_pos
  return PMT_POSITION


def xyz_to_phi_theta(x, y, z):
   phi = math.atan2(y, x)
   r = (x**2 + y**2 + z**2)**.5
   theta = math.acos(z / r)
   return phi, theta

#change directory
class cd:
   '''
   Context manager for changing the current working directory
   '''
   def __init__(self, newPath):
      self.newPath = newPath

   def __enter__(self):
      self.savedPath = os.getcwd()
      os.chdir(self.newPath)

   def __exit__(self, etype, value, traceback):
      os.chdir(self.savedPath)

def tof(pmt_position, vertex_position):

  # ceff = 16.95 #cm/ns From KatLTVertex.cc source code in Kat

  return np.linalg.norm(pmt_position-vertex_position,2) * 100.0/16.95




# Convert the phi theta information to row and column index in 2D grid
def phi_theta_to_row_col(phi, theta, rows=ROWS, cols=COLS):
   # phi is in [-pi, pi], theta is in [0, pi]
   row = min(rows/2 + (math.floor((rows/2)*phi/math.pi)), rows-1)
   row = max(row, 0)
   col = min(math.floor(cols*theta/math.pi), cols-1);
   col = max(col, 0)
   return int(row), int(col)

# Calculating the angle between two input vectors
def calculate_angle(vec1, vec2):
  x1,y1,z1 = vec1
  x2,y2,z2 = vec2
  inner_product = x1*x2 + y1*y2 + z1*z2
  len1 = (x1**2 + y1**2 + z1**2)**0.5
  len2 = (x2**2 + y2**2 + z2**2)**0.5
  return math.acos(float(inner_product)/float(len1*len2))

# Converting Cartesian position to 2D Grid
def xyz_to_row_col(pmt_index, PMT_POSITION,rows=ROWS, cols=COLS):
   x, y, z = tuple(PMT_POSITION[pmt_index])
   return phi_theta_to_row_col(*xyz_to_phi_theta(x, y, z), rows=rows, cols=cols)


# Set up the clocl to start ticking on the first incoming photon of a given events.
def set_clock(tree, evt):
  tree.GetEntry(evt)
  time_array = []
  for i in range(tree.N_phot):
    time_array.append(tree.PE_time[i])
  return clock(np.array(time_array).min())

# Save input file as a .json file.
def savefile(saved_file, appendix, filename, pathname):
    if not os.path.exists(pathname):
     os.mkdir(pathname)
    with cd(pathname):
        with open(filename, 'w') as datafile:
          json.dump(saved_file, datafile)

def calculate_tzero(t, tof, charge):
  return np.sum((t - tof) * charge)/ np.sum(charge)

def append_file(key, val, filename_array):
    if key not in filename_array.keys():
        filename_array[key] = [val]
    else:
        filename_array[key].append(val)
    return filename_array

def get_name(file):
     '''
     Retrieve run number from file name
     '''
     zdab_regex = re.compile(r"^eventfile_sph_out_(.+)_[a-z].*_1k_(\d+)\.\d+\.\d+.pickle$")
     matches = zdab_regex.match(file)
     if matches:
            return str(matches.group(1)), str(matches.group(2))
     else:
            return 0