#154: Merged with master

GetLLM/GetLLM.py

@@ -189,10 +189,10 @@ def _parse_args():
                     default=USE_ONLY_THREE_BPMS_FOR_BETA_FROM_PHASE, dest="use_only_three_bpms_for_beta_from_phase")
     parser.add_option("-j", "--numbpm",
                     help="Number of different BPM combinations for beta-calculation, default = 10",
-                    metavar="NUMBER_OF_BPMS", default=NUMBER_OF_BPMS, dest="number_of_bpms")
+                    metavar="NUMBER_OF_BPMS", type=int, default=NUMBER_OF_BPMS, dest="number_of_bpms")
     parser.add_option("-i", "--range",
                     help="Range of BPM for beta-calculation (>=3 and odd), default = 11",
-                    metavar="RANGE_OF_BPMS", default=RANGE_OF_BPMS, dest="range_of_bpms")
+                    metavar="RANGE_OF_BPMS", type=int, default=RANGE_OF_BPMS, dest="range_of_bpms")
     parser.add_option("-r", "--average_tune",
                     help="Set to 1 to use average tune for all BPMs instead of specific for each one.",
                     metavar="AVERAGE_TUNE", type="int", default=AVERAGE_TUNE, dest="use_average")
@@ -791,7 +791,7 @@ def _analyse_src_files(getllm_d, twiss_d, files_to_analyse, nonlinear, turn_by_t
             if use_average:
                 twiss_file_x.MUX = twiss_file_x.AVG_MUX
             if calibration_twiss is not None:
-                twiss_file_x.AMPX, twiss_file_x.ERRAMPX = _get_calibrated_amplitudes(twiss_file_x, calibration_twiss, "X")
+                twiss_file_x.AMPX, twiss_file_x.ERRAMPX, twiss_file_x.CALIBRATION, twiss_file_x.ERROR_CALIBRATION = _get_calibrated_amplitudes(twiss_file_x, calibration_twiss, "X")
             try:
                 dppi = getattr(twiss_file_x, "DPP", 0.0)
             except AttributeError:
@@ -869,7 +869,7 @@ def _analyse_src_files(getllm_d, twiss_d, files_to_analyse, nonlinear, turn_by_t
             if use_average:
                 twiss_file_y.MUY = twiss_file_y.AVG_MUY
             if calibration_twiss is not None:
-                twiss_file_y.AMPY, twiss_file_y.ERRAMPY = _get_calibrated_amplitudes(twiss_file_y, calibration_twiss, "Y")
+                twiss_file_y.AMPY, twiss_file_y.ERRAMPY, twiss_file_y.CALIBRATION, twiss_file_y.ERROR_CALIBRATION = _get_calibrated_amplitudes(twiss_file_y, calibration_twiss, "Y")
             try:
                 dppi = getattr(twiss_file_y, "DPP", 0.0)
             except AttributeError:
@@ -1155,20 +1155,47 @@ def _calculate_kick(kick_times, getllm_d, twiss_d, phase_d, beta_d, mad_twiss, m
 
 
 def _get_calibrated_amplitudes(drive_file, calibration_twiss, plane):
-    calibration_file = calibration_twiss[plane]
-    cal_amplitudes = []
-    err_cal_amplitudes = []
-    for bpm_name in drive_file.NAME:
-        drive_index = drive_file.indx[bpm_name]
-        cal_amplitude = getattr(drive_file, "AMP" + plane)[drive_index]
-        err_cal_amplitude = 0.
-        if bpm_name in calibration_file.NAME:
-            cal_index = calibration_file.indx[bpm_name]
-            cal_amplitude = cal_amplitude * calibration_file.CALIBRATION[cal_index]
-            err_cal_amplitude = calibration_file.ERROR_CALIBRATION[cal_index]
-        cal_amplitudes.append(cal_amplitude)
-        err_cal_amplitudes.append(err_cal_amplitude)
-    return array(cal_amplitudes), array(err_cal_amplitudes)
+   calibration_file = calibration_twiss[plane]
+   cal_amplitudes = []
+   err_cal_amplitudes = []
+   calibration_value = []
+   calibration_error = []
+   if plane == "X":
+      tune = "Q1"
+      tune_rms = "Q1RMS"
+      natural_tune = "NATQ1"
+      natural_tune_rms = "NATQ1RMS"
+   elif plane == "Y":
+      tune = "Q2"
+      tune_rms = "Q2RMS"
+      natural_tune = "NATQ2"
+      natural_tune_rms = "NATQ2RMS"
+   tune_value = getattr(drive_file,tune)
+   #print "DRIVE TUNE"
+   #print tune_value
+   #print tune_value_rms
+   #print "NATURAL TUNE"
+   #print natural_tune_value
+   #print natural_tune_value_rms
+   tune_value_rms = getattr(drive_file,tune_rms)
+   natural_tune_value = getattr(drive_file,natural_tune_rms)
+   natural_tune_value_rms = getattr(drive_file,natural_tune_rms)
+   for bpm_name in drive_file.NAME:
+       drive_index = drive_file.indx[bpm_name]
+       cal_amplitude = getattr(drive_file, "AMP" + plane)[drive_index]
+       err_cal_amplitude = 0.
+       if bpm_name in calibration_file.NAME:
+           cal_index = calibration_file.indx[bpm_name]
+           cal_amplitude = cal_amplitude * calibration_file.CALIBRATION[cal_index]
+           err_cal_amplitude = cal_amplitude * calibration_file.ERROR_CALIBRATION[cal_index]
+           calibration_value.append(calibration_file.CALIBRATION[cal_index])
+           calibration_error.append(calibration_file.ERROR_CALIBRATION[cal_index])
+       else:
+           calibration_value.append(1.)
+           calibration_error.append(0.)
+       cal_amplitudes.append(cal_amplitude)
+       err_cal_amplitudes.append(err_cal_amplitude)
+   return array(cal_amplitudes), array(err_cal_amplitudes),array(calibration_value),array(calibration_error)
 # END _get_calibrated_amplitudes --------------------------------------------------------------------
 
 

GetLLM/algorithms/compensate_ac_effect.py

@@ -360,7 +360,7 @@ def get_free_phase_eq(MADTwiss, Files, Qd, Q, psid_ac2bpmac, plane, bd, op, Qmdl
     return result, muave, bpm
 
 
-def get_free_beta_from_amp_eq(MADTwiss_ac, Files, Qd, Q, psid_ac2bpmac, plane, bd, op):
+def get_free_beta_from_amp_eq(MADTwiss_ac, Files, Qd, Q, psid_ac2bpmac, plane, bd, op):#,calibration,calibration_error):
     #-- Select common BPMs
     all_bpms = utils.bpm.model_intersect(
         utils.bpm.intersect(Files),
@@ -373,7 +373,7 @@ def get_free_beta_from_amp_eq(MADTwiss_ac, Files, Qd, Q, psid_ac2bpmac, plane, b
     print ("skowron: Please fix me !!!! ")
     print ("skowron: op is sometimes the machine name and sometimes lhcphase1 flag  ")
     print "op"
-    print op
+    print op 
     if op == "1": 
         print "here"
         good_bpms_for_kick = intersect_bpm_list_with_arc_bpms(bpms)
@@ -424,11 +424,28 @@ def get_free_beta_from_amp_eq(MADTwiss_ac, Files, Qd, Q, psid_ac2bpmac, plane, b
 
     #-- Loop for files
     betall = np.zeros((len(all_bpms), len(Files)))
+    betall_err = np.zeros((len(all_bpms), len(Files)))
     for i in range(len(Files)):
         if plane == 'H':
             amp = np.array(
                 [2 * Files[i].AMPX[Files[i].indx[b[1]]] for b in all_bpms]
             )
+            try:
+                 amp_err = np.array(
+                    [Files[i].ERRAMPX[Files[i].indx[b[1]]] for b in all_bpms]
+                )
+                 print "CALIBRATION"
+                 print Files[i].CALIBRATION
+                 calibration = np.array(
+                    [Files[i].CALIBRATION[Files[i].indx[b[1]]] for b in all_bpms]
+                )
+                 calibration_error = np.array(
+                    [Files[i].ERROR_CALIBRATION[Files[i].indx[b[1]]] for b in all_bpms]
+                )
+            except AttributeError:
+                  amp_err = np.zeros(len(amp))
+                  calibration = np.ones(len(amp))
+                  calibration_error = np.zeros(len(amp))
             psid = bd * 2 * np.pi * np.array(
                 [Files[i].MUX[Files[i].indx[b[1]]] for b in all_bpms]
             )  # bd flips B2 phase to B1 direction
@@ -439,8 +456,22 @@ def get_free_beta_from_amp_eq(MADTwiss_ac, Files, Qd, Q, psid_ac2bpmac, plane, b
             psid = bd * 2 * np.pi * np.array(
                 [Files[i].MUY[Files[i].indx[b[1]]] for b in all_bpms]
             )  # bd flips B2 phase to B1 direction
-
-        # This loop is just to fix the phase jump at the beginning of the ring.
+            try:
+                 amp_err = np.array(
+                    [Files[i].ERRAMPY[Files[i].indx[b[1]]] for b in all_bpms]
+                )
+                 calibration = np.array(
+                    [Files[i].CALIBRATION[Files[i].indx[b[1]]] for b in all_bpms]
+                )
+                 calibration_error = np.array(
+                    [Files[i].ERROR_CALIBRATION[Files[i].indx[b[1]]] for b in all_bpms]
+                )
+            except AttributeError:
+                  amp_err = np.zeros(len(amp))
+                  calibration = np.ones(len(amp))
+                  calibration_error = np.zeros(len(amp))
+
+# This loop is just to fix the phase jump at the beginning of the ring.
         for k in range(len(all_bpms)):
             try:
                 if all_bpms[k][0] > s_lastbpm:
@@ -453,17 +484,29 @@ def get_free_beta_from_amp_eq(MADTwiss_ac, Files, Qd, Q, psid_ac2bpmac, plane, b
         Psid[k_bpmac:] = Psid[k_bpmac:] - 2 * np.pi * Qd
         bet = ((amp / sqrt2j[i]) ** 2 *
                (1 + r ** 2 + 2 * r * np.cos(2 * Psid)) / (1 - r ** 2))
+        print "BETA ERROR CALIBRATION"
+        print amp_err
+        #print bet_err
+        #bet_err = (2 * calibration * bet * calibration_error)**2 
+        bet_err = ((amp_err / sqrt2j[i]) ** 2 *
+               (1 + r ** 2 + 2 * r * np.cos(2 * Psid)) / (1 - r ** 2))
+        print "BETA ERROR"
+        print bet_err
+        print bet_err
         for bpm_index in range(len(all_bpms)):
             betall[bpm_index][i] = bet[bpm_index]
-
+            # betall_err[bpm_index][i] = bet_err[bpm_index]
+            betall_err[bpm_index][i] = 2 / calibration[bpm_index] * bet[bpm_index] * calibration_error[bpm_index]
     #-- Output
     result = {}
     bb = []
     for k in range(len(all_bpms)):
         betave = np.mean(betall[k])
+        bet_errave = np.mean(betall_err[k])
         betstd = np.std(betall[k])
+        bet_err_total = (bet_errave**2+betstd**2)**0.5
         bb.append((betave - betmdl[k]) / betmdl[k])
-        result[all_bpms[k][1]] = [betave, betstd, all_bpms[k][0]]
+        result[all_bpms[k][1]] = [betave, bet_err_total, all_bpms[k][0]]
     bb = math.sqrt(np.mean(np.array(bb) ** 2))
 
     return result, bb, all_bpms

GetLLM/get_psb_calibration_factors.py

@@ -0,0 +1,81 @@
+import sys
+import os
+from os.path import join, dirname, abspath, pardir
+from optparse import OptionParser
+from collections import OrderedDict
+import numpy as np
+import pandas as pd
+
+new_path = abspath(join(dirname(abspath(__file__)), pardir))
+if new_path not in sys.path:
+    sys.path.append(new_path)
+
+import tfs_files.tfs_file_writer as tfs_writer
+from tfs_files import tfs_pandas
+
+
+def _parse_args():
+    parser = OptionParser()
+    parser.add_option("-r", "--results",
+                    help="results folder",
+                    metavar="results",dest="results")
+    parser.add_option("-o", "--output_path",
+                    help="output path",
+                    metavar="OUTPUT",dest="output_path")
+    options, _ = parser.parse_args()
+    return options.results,options.output_path
+
+
+
+def get_beta_from_amplitude(beta_amplitude,plane):
+    beta_file_amplitude = tfs_pandas.read_tfs(os.path.join(beta_amplitude,"getampbeta" + str(plane.lower()) + "_free.out"))
+    return beta_file_amplitude
+
+
+def get_beta_from_phase(beta_phase,plane):
+    beta_file_phase = tfs_pandas.read_tfs(os.path.join(beta_phase,"getbeta" + str(plane.lower()) + "_free.out"))
+    return beta_file_phase
+
+
+def get_betas_values_with_errors(beta_file_phase,beta_file_amplitude,plane):
+    beta_error_phase = "ERRBET" + str(plane.upper())
+    beta_error_amplitude = "BET" + str(plane.upper()) + "STD"
+    beta_model = "BET" + str(plane.upper()) + "MDL"
+    beta = "BET" + str(plane.upper())
+    beta_phase = getattr(beta_file_phase,beta)
+    beta_model = getattr(beta_file_phase,beta_model)
+    beta_amplitude = getattr(beta_file_amplitude,beta)
+    beta_error_phase_values = getattr(beta_file_phase,beta_error_phase)
+    beta_error_amplitude_values = getattr(beta_file_amplitude,beta_error_amplitude)
+    names = beta_file_phase.NAME
+    position = beta_file_phase.S
+    return beta_phase,beta_error_phase_values,beta_amplitude,beta_error_amplitude_values,beta_model,names,position
+
+
+def get_calibration_factors(beta_phase,beta_amplitude,beta_error_phase,beta_error_amplitude):
+    beta_calibration = (beta_phase/beta_amplitude)**0.5
+    beta_calibration_error = (0.25*(beta_phase/beta_amplitude**3)*beta_error_amplitude**2 + 0.25*(1/(beta_phase*beta_amplitude))*beta_error_phase**2)**0.5
+    return beta_calibration,beta_calibration_error
+
+
+def writting_calibration_factors_to_file(name,position,beta_calibration,beta_calibration_error,file_output):
+    tfs_calibration_pandas = file_output
+    all_summary = pd.DataFrame(OrderedDict([('NAME', name),('S', position),('CALIBRATION',beta_calibration),('ERROR_CALIBRATION',beta_calibration_error)]))
+    tfs_pandas.write_tfs(tfs_calibration_pandas,all_summary,save_index=False)
+
+
+def main(results_directory,output_directory):
+    planes = ["x","y"] 
+    for plane in planes:
+        file_output_name = "calibration_" + str(plane.lower()) + ".out"
+        file_output_results = os.path.join(output_directory,file_output_name)
+        beta_from_phase_file = get_beta_from_phase(results_directory,str(plane.lower())) 
+        beta_from_amplitude_file = get_beta_from_amplitude(results_directory,str(plane.lower()))
+        [beta_phase, beta_error_phase, beta_amplitude, beta_error_amplitude, beta_model, name, position] = get_betas_values_with_errors(beta_from_phase_file,beta_from_amplitude_file,plane)
+        [beta_calibration, beta_calibration_error] = get_calibration_factors(beta_phase,beta_amplitude,beta_error_phase,beta_error_amplitude)
+        writting_calibration_factors_to_file(name, position, beta_calibration, beta_calibration_error, file_output_results)
+
+
+if __name__ == "__main__":
+    (_file_results,_output_path) = _parse_args()
+    main(_file_results,_output_path)

GetLLM/getsuper.py

@@ -662,11 +662,13 @@ def _get_tunes(model_file, fileslist):
 
 def linreg(X, Y):
     """
-    Summary
+    Returns coefficients to the regression line "y=ax+b" from x[] and y[], and R^2 Value
+
+    Summary:
         Linear regression of y = ax + b
-    Usage
+
+    Usage:
         real, real, real = linreg(list, list)
-    Returns coefficients to the regression line "y=ax+b" from x[] and y[], and R^2 Value
     """
     if len(X) != len(Y):  raise ValueError, 'unequal length'
     N = len(X)