From cdb6cc3c1421b514631d3410520f9cf33b1a7ddd Mon Sep 17 00:00:00 2001
From: Ryoichi Miyamoto <ryoichi.miyamoto@esss.se>
Date: Fri, 29 Jun 2018 17:33:40 +0200
Subject: [PATCH] A stand-alone script to convert the IBSimu's source output to
the TraceWin's format.
---
ess/ibsimu2tw.py | 122 +++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 122 insertions(+)
create mode 100644 ess/ibsimu2tw.py
diff --git a/ess/ibsimu2tw.py b/ess/ibsimu2tw.py
new file mode 100644
index 0000000..9e6c108
--- /dev/null
+++ b/ess/ibsimu2tw.py
@@ -0,0 +1,122 @@
+
+'''
+ - Convert an IBSimu's output (export_path_manager_data) to a TraceWin's binary.
+ - Use the ipynb version for checking.
+ - IBSimu's output format:
+
+ 0: index
+ 1: x [m]
+ 2: x' [rad]
+ 3: y [m]
+ 4: y' [rad]
+ 5: phi [rad]
+ 6: (pz-p0)/p0 [1]
+ 9: mc2 [GeV]
+
+ R. Miyamoto (2018.06.25)
+'''
+
+
+#-- Constants and inputs
+
+mc2_H1 = 938.272
+freq = 352.21
+
+file_name_in = 'tracewin_path_format_2d_10M.txt'
+file_name_out_H1 = 'ISRC.H1_ibsimu.2d_10M.83mA.dst'
+file_name_out_H2 = 'ISRC.H2_ibsimu.2d_10M.83mA.dst' # Comment out if H2 file isn't needed
+file_name_out_H3 = 'ISRC.H3_ibsimu.2d_10M.83mA.dst' # Comment out if H3 file isn't needed
+
+#W_cut = 0.01 # Filter for W [MeV]. Only for experts.
+
+##################################################
+# No need to change below
+##################################################
+
+#-- Libraries
+
+import numpy as np
+import sys
+
+import matplotlib.pyplot as plt
+
+#-- Function for saving
+
+def save_dst(file_name, x6, Ibeam, mc2=938.272, freq=352.21):
+ '''
+ - x6 dimensions can be either (Nptcl, 6) or (Nptcl*6)
+ - Coordinates: ( x, x', y, y', phi, W)
+ - Units : ( cm, rad, cm, rad, rad, MeV)
+ '''
+ import numpy, struct
+
+ x6 = numpy.reshape(x6, (numpy.size(x6)))
+ with open(file_name, 'wb') as file:
+ out = struct.pack('b' , 125 )
+ out += struct.pack('b' , 100 )
+ out += struct.pack('i' , len(x6))
+ out += struct.pack('d' , Ibeam )
+ out += struct.pack('d' , freq )
+ out += struct.pack('b' , 125 ) # Typo in the manual !!!!
+ out += struct.pack('%s'%len(x6)+'d', *x6 )
+ out += struct.pack('d' , mc2 )
+ file.write(out)
+
+#-- Parameters from the header
+
+with open(file_name_in, 'r') as file:
+ for i, lin in enumerate(file):
+ if 'beam' in lin: Ibeam = float(lin.partition(':')[-1].split()[0])*1e3 # A => mA
+ if 'MOMENTUM' in lin: cp = float(lin.split()[0] )*1e3 # GeV => MeV
+ if i == 8 : break
+
+#-- Read data (could take tens of seconds)
+
+data = np.loadtxt(file_name_in, skiprows=8)
+
+#-- Adjust coordinates and units
+
+data[:,1] *= 1e2 # x : m => cm
+data[:,3] *= 1e2 # y : m => cm
+data[:,9] *= 1e3 # mc2: GeV => MeV
+data[:,6] = np.sqrt(data[:,9]**2+(1.0+data[:,2]**2+data[:,4]**2)*(cp*(1.0+data[:,6]))**2)-data[:,9] # dpz => W
+
+try : data = data[np.where(np.abs(data[:,6]-0.075) < W_cut)]
+except: pass
+
+#-- Separate species
+
+mc2_lst = sorted(list(set(list(data[:,9]))))
+
+x6_H1 = data[np.where(data[:,9] == mc2_lst[0])]
+x6_H1 = x6_H1[:,1:7]
+Ibeam_H1 = Ibeam*len(x6_H1)/len(data)
+
+try:
+ x6_H2 = data[np.where(data[:,9] == mc2_lst[1])]
+ x6_H2 = x6_H2[:,1:7]
+ Ibeam_H2 = Ibeam*len(x6_H2)/len(data)
+except:
+ pass
+
+try:
+ x6_H3 = data[np.where(data[:,9] == mc2_lst[2])]
+ x6_H3 = x6_H3[:,1:7]
+ Ibeam_H3 = Ibeam*len(x6_H3)/len(data)
+except:
+ pass
+
+#-- Save outputs
+
+save_dst( file_name_out_H1, x6_H1, Ibeam_H1, 1.0*mc2_H1, freq)
+
+try : save_dst(file_name_out_H2, x6_H2, Ibeam_H2, 2.0*mc2_H1, freq)
+except: pass
+try : save_dst(file_name_out_H3, x6_H3, Ibeam_H3, 3.0*mc2_H1, freq)
+except: pass
+
+#-- Ending
+
+sys.exit()
+
+#--
--
GitLab