corrections according to the reviewer's comments

examples/sim/efield2voltage_manually.py

@@ -9,7 +9,7 @@
 import numpy as np
 from grand.sim.efield2voltage import Efield2Voltage
 from grand.sim.noise.galaxy import galactic_noise
-import grand.sim.detector.rf_chain2 as grfc
+import grand.sim.detector.rf_chain as grfc
 
 rng = np.random.default_rng(0)     
 

grand/geo/coordinates.py

@@ -102,22 +102,6 @@ def geoid_undulation(latitude=None, longitude=None):
 # Define functions to transform from one coordinate representation to
 # another coordinate representation. Cartesian, Spherical, and Horizontal
 # coordinate representation are defined.
-
-#The old angular convention
-#def _cartesian_to_spherical(
-#    x: Union[float, int, np.ndarray],
-#    y: Union[float, int, np.ndarray],
-#    z: Union[float, int, np.ndarray],
-#) -> Union[Tuple[float, ...], Tuple[np.ndarray, ...]]:
-#    """Transform Cartesian coordinates to spherical"""
-#    rho2 = x ** 2 + y ** 2
-#    rho = np.sqrt(rho2)
-#    theta = np.rad2deg(np.arctan2(rho, z))
-#    phi = np.rad2deg(np.arctan2(y, x))
-#    r = np.sqrt(rho2 + z ** 2)
-#
-#    return theta, phi, r
-
 def _cartesian_to_spherical(
     x: Union[float, int, np.ndarray],
     y: Union[float, int, np.ndarray],
@@ -126,16 +110,11 @@ def _cartesian_to_spherical(
     """Transform Cartesian coordinates to spherical"""
     rho2 = x ** 2 + y ** 2
     rho = np.sqrt(rho2)
-    theta = 180. - np.rad2deg(np.arctan2(rho, z))
-    phi = np.rad2deg(np.arctan2(y, x)) + 180.
+    theta = np.rad2deg(np.arctan2(rho, z))
+    phi = np.rad2deg(np.arctan2(y, x))
     r = np.sqrt(rho2 + z ** 2)
-    if phi==360:
-        phi=0
-    else:
-        phi=phi
-    return theta, phi, r
-
 
+    return theta, phi, r
 
 
 # Horizontal has an axis fixed to geographic North, so it can not be
@@ -150,38 +129,22 @@ def _cartesian_to_horizontal(
     theta, phi, r = _cartesian_to_spherical(x, y, z)
     return _spherical_to_horizontal(theta, phi, r)
 
-#The old angular convention
-#def _spherical_to_cartesian(
-#    theta: Union[float, int, np.ndarray],
-#    phi: Union[float, int, np.ndarray],
-#    r: Union[float, int, np.ndarray],
-#) -> Union[Tuple[float, ...], Tuple[np.ndarray, ...]]:
-#    """Transform spherical coordinates to Cartesian"""
-#    cos_theta = np.cos(np.deg2rad(theta))
-#    sin_theta = np.sin(np.deg2rad(theta))
-#
-#    x = r * np.cos(np.deg2rad(phi)) * sin_theta
-#    y = r * np.sin(np.deg2rad(phi)) * sin_theta
-#    z = r * cos_theta
-#
-#    return x, y, z
 
 def _spherical_to_cartesian(
     theta: Union[float, int, np.ndarray],
     phi: Union[float, int, np.ndarray],
     r: Union[float, int, np.ndarray],
 ) -> Union[Tuple[float, ...], Tuple[np.ndarray, ...]]:
     """Transform spherical coordinates to Cartesian"""
-    cos_theta = np.cos(np.deg2rad(180. - theta))
-    sin_theta = np.sin(np.deg2rad(180. - theta))
+    cos_theta = np.cos(np.deg2rad(theta))
+    sin_theta = np.sin(np.deg2rad(theta))
 
-    x = r * np.cos(np.deg2rad(phi + 180.)) * sin_theta
-    y = r * np.sin(np.deg2rad(phi + 180.)) * sin_theta
+    x = r * np.cos(np.deg2rad(phi)) * sin_theta
+    y = r * np.sin(np.deg2rad(phi)) * sin_theta
     z = r * cos_theta
 
     return x, y, z
-
-
+
 
 def _spherical_to_horizontal(
     theta: Union[float, int, np.ndarray],
@@ -190,8 +153,8 @@ def _spherical_to_horizontal(
 ) -> Tuple[Union[float, np.ndarray], Union[float, np.ndarray], Union[float, np.ndarray]]:
     """Transform spherical coordinates to horizontal"""
     # return 0.5 * np.pi - phi, 0.5 * np.pi - theta, r
-    #return 90.0 - phi, 90.0 - theta, r #(the old angular convention)
-    return 270.0 - phi, theta - 90.0, r
+    return 90.0 - phi, 90.0 - theta, r
+
 
 # Horizontal has an axis fixed to geographic North, so it can not be
 # converted like cartesian and spherical. If conversion is done, then
@@ -213,8 +176,8 @@ def _horizontal_to_spherical(
 ) -> Tuple[Union[float, np.ndarray], Union[float, np.ndarray], Union[float, np.ndarray]]:
     """Transform horizontal coordinates to spherical"""
     # return 0.5 * np.pi - elevation, 0.5 * np.pi - azimuth, norm
-    #return 90.0 - elevation, 90.0 - azimuth, norm #(the old angular convention)
-    return 90.0 + elevation, 270.0 - azimuth, norm
+    return 90.0 - elevation, 90.0 - azimuth, norm
+
 
 # -----------Base Representation------------
 class Coordinates(np.ndarray):
@@ -1297,7 +1260,6 @@ def __init__(
             height=height,  # height of LTP's location/origin
             location=location,  # location of LTP in Geodetic, GRANDCS, or ECEF
             orientation="NWU",  # orientation of LTP. 'NWU', 'ENU' etc
-            #orientation="SED",  # orientation of LTP. 'NWU', 'ENU' etc
             magnetic=True,  # shift orientation by magnetic declination?
             magmodel="IGRF13",  # if shift, which magnetic model to use?
             declination=None,  # or simply provide the magnetic declination
@@ -1321,4 +1283,4 @@ def grandcs_to_ltp(self, ltp):
 
 # RK TODO: Rework on this class.
 class Rotation(_Rotation):
-    pass
+    pass

grand/sim/detector/antenna_model.py

@@ -1,3 +1,11 @@
+# Modified by SN in April 2024, in order to include the models for the effective length of the antennas produced by the codes NEC and matlab from HOU group. 
+# The coresponding files are Light_GP300Antenna_nec_*_leff.npz and Light_GP300Antenna_mat_*_leff.npz for NEC and matlab respectively (* stands for the 3 arms X-NS, Y-EW and Z-vertical)
+# The default effective length model were produced using HFSS simulation code by Xidian group. 
+# The coresponding files of the default model are Light_GP300Antenna_EWarm_leff.npz, Light_GP300Antenna_SNarm_leff.npz and Light_GP300Antenna_Zarm_leff.npz.
+# For details on how these simulations were produced using the different codes and comparisons between different model (HFSS, NEC matlab) please refer to the documents and presentations 
+# posted in forge software forum https://forge.in2p3.fr/projects/software-forum/dmsf (for NEC and matlab) and https://forge.in2p3.fr/documents/1353, https://forge.in2p3.fr/documents/1354 for the HFSS simulations.
+
+
 from logging import getLogger
 
 from grand import grand_add_path_data

scripts/T1_trigger_offline.py

@@ -0,0 +1,135 @@
+#!/pbs/home/x/xtian/.conda/envs/grandlib2304/bin/python3.9
+import numpy as np
+import sys
+import grand.dataio.root_trees as rt
+
+def extract_trigger_parameters(trace, trigger_config, baseline=0):
+    # Extract the trigger infos from a trace
+
+    # Parameters :
+    # ------------
+    # trace, numpy.ndarray: 
+    # traces in ADC unit
+    # trigger_config, dict:
+    # the trigger parameters set in DAQ
+
+    # Returns :
+    # ---------
+    # Index in the trace when the first T1 crossing happens
+    # Indices in the trace of T2 crossing happens
+    # Number of T2 crossings
+    # Q, Peak/NC
+
+    # Find the position of the first T1 crossing
+    index_t1_crossing = np.where((trace) > trigger_config["th1"],
+                                 np.arange(len(trace)), -1)
+    dict_trigger_infos = dict()
+    mask_T1_crossing = (index_t1_crossing != -1)
+    if sum(mask_T1_crossing) == 0:
+        # No T1 crossing 
+        raise ValueError("No T1 crossing!")
+    dict_trigger_infos['index_T1_crossing'] = None
+    # Tquiet to decide the quiet time before the T1 crossing 
+    for i in index_t1_crossing[mask_T1_crossing]:
+       # Abs value not exceeds the T1 threshold
+       if i - trigger_config["t_quiet"]//2 < 0:
+          raise ValueError("Not enough data before T1 crossing!")
+       if np.all((trace[np.max([0, i - trigger_config['t_quiet'] // 2]):i]) <= trigger_config["th1"]):
+          dict_trigger_infos["index_T1_crossing"] = i
+          # the first T1 crossing satisfying the quiet condition
+          break
+    if dict_trigger_infos['index_T1_crossing'] == None:
+       raise ValueError("No T1 crossing with Tquiet satified!")
+    # The trigger logic works for the timewindow given by T_period after T1 crossing.
+    # Count number of T2 crossings, relevant pars: T2, NCmin, NCmax, T_sepmax
+    # From ns to index, divided by two for 500MHz sampling rate
+    period_after_T1_crossing = trace[dict_trigger_infos["index_T1_crossing"]:dict_trigger_infos["index_T1_crossing"]+trigger_config['t_period']//2]
+    # All the points above +T2
+    positive_T2_crossing = (np.array(period_after_T1_crossing) > trigger_config['th2']).astype(int)
+    # Positive crossing, the point before which is below T2.
+    mask_T2_crossing_positive = np.diff(positive_T2_crossing) == 1
+    # if np.sum(mask_T2_crossing_positive) > 0:
+    #     index_T2_crossing_positive = np.arange(len(period_after_T1_crossing) - 1)[mask_T2_crossing_positive]
+    negative_T2_crossing = (np.array(period_after_T1_crossing) < - trigger_config['th2']).astype(int)
+    mask_T2_crossing_negative = np.diff(negative_T2_crossing) == 1
+    # if np.sum(mask_T2_crossing_negative) > 0:
+    #     index_T2_crossing_negative = np.arange(len(period_after_T1_crossing) - 1)[mask_T2_crossing_negative]
+    # n_T2_crossing_negative = np.len(index_T2_crossing_positive)
+    # Register the first T1 crossing as a T2 crossing
+    mask_first_T1_crossing = np.zeros(len(period_after_T1_crossing), dtype=bool)
+    mask_first_T1_crossing[0] = True
+    # mask_first_T1_crossing[1:] = (mask_T2_crossing_positive | mask_T2_crossing_negative)
+    mask_first_T1_crossing[1:] = (mask_T2_crossing_positive)
+    index_T2_crossing = np.arange(len(period_after_T1_crossing))[mask_first_T1_crossing]
+    n_T2_crossing = 1 # Starting from the first T1 crossing.
+    dict_trigger_infos["index_T2_crossing"] = [0]
+    if len(index_T2_crossing) > 1:
+      for i, j in zip(index_T2_crossing[:-1], index_T2_crossing[1:]):
+          # The separation between successive T2 crossings
+          time_separation = (j - i) * 2
+          if time_separation < trigger_config["t_sepmax"]:
+              n_T2_crossing += 1
+              dict_trigger_infos["index_T2_crossing"].append(j)
+          else:
+              # Violate the maximum separation, fail to trigger
+              raise ValueError(f"Violating Tsepmax, the separation is {time_separation} ns.")
+    else:
+      n_T2_crossing = 1
+      j = 1
+    # Change the reference of indices of T2 crossing
+    dict_trigger_infos["index_T2_crossing"] = np.array(dict_trigger_infos["index_T2_crossing"]) + dict_trigger_infos["index_T1_crossing"]
+    dict_trigger_infos["NC"] = n_T2_crossing
+    # Calulate the peak value
+    dict_trigger_infos["Q"] = (np.max(np.abs(period_after_T1_crossing[:j])) - baseline) / dict_trigger_infos["NC"]
+    return dict_trigger_infos
+
+dict_trigger_parameter = dict([
+  ("t_quiet", 512),
+  ("t_period", 512),
+  ("t_sepmax", 10),
+  ("nc_min", 2),
+  ("nc_max", 8),
+  ("q_min", 0),
+  ("q_max", 255),
+  ("th1", 100),
+  ("th2", 50),
+  # Configs of readout timewindow
+  ("t_pretrig", 960),
+  ("t_overlap", 64),
+  ("t_posttrig", 1024)
+  ])
+
+
+if __name__ == "__main__":
+  # Read the traces from experimental data
+  fname = sys.argv[1]
+  file = rt.DataFile(fname)
+  n_entries = file.tadc.get_number_of_entries()
+  # Pad zeros at the head of the trace to statify the Tquiet condition
+  zero_head = np.zeros(dict_trigger_parameter["t_quiet"] // 2, dtype=int)
+
+  trigger_index = []
+  # Loop over all entries
+  for k in range(n_entries):
+    file.tadc.get_entry(k)
+    # Loop over four channels
+    #for v in range(4):
+    for v in range(3):
+      trace = file.tadc.trace_ch[0][v]
+      # Zero padding at first
+      # trace_padded = np.concatenate((zero_head, trace))
+      try:
+        # Check if trigger
+        trigger_infos = extract_trigger_parameters(trace, dict_trigger_parameter)
+        # Save the triggered traces
+        if trigger_infos["NC"] >= dict_trigger_parameter["nc_min"] and trigger_infos["NC"] <= dict_trigger_parameter["nc_max"]:
+          trigger_index.append(k)
+          break
+      except ValueError:
+        # No T1 crossing, no trigger
+        # print(k, ": No trigger.")
+        pass
+
+  print(f"{fname}: {len(trigger_index)} out of {n_entries} triggered.")
+  if len(trigger_index) > 0:
+    np.savetxt(f"./{fname.split('/')[-1]}.trigger.txt", trigger_index, delimiter=', ', fmt='%d', header=str(dict_trigger_parameter))

sim2root/CoREASRawRoot/CoreasToRawROOT.py

@@ -178,11 +178,9 @@ def CoreasToRawRoot(file, simID=None):
   MesonEnergyCut    = HadronEnergyCut # mesons are hadronic, so this should be fine
 
   #parallel = read_list_of_params(inp_input, "PARALLEL") # COREAS-only
-  #ECTCUT = parallel[0]
-  #ECTMAX = parallel[1]
-  ECTCUT = 1.
-  ECTMAX = 1.
-
+  ECTCUT = parallel[0]
+  ECTMAX = parallel[1]
+
   # PARALLEL = [ECTCUT, ECTMAX, MPIID, FECTOUT]
   # ECTCUT: limit for subshowers GeV
   # ECTMAX: maximum energy for complete shower GeV

sim2root/CoREASRawRoot/CorsikaInfoFuncs.py

@@ -227,7 +227,7 @@ def read_long(pathLongFile):
 
         for line in file:
             # use a regex to search for a minus sign that is not part of an exponent
-            if search(r"(?<!E)(-)(?=\d)", line):
+            if search(r"(?<!e)(-)(?=\d)", line):
                 # if the minus sign is not part of an exponent, replace it with a space and a minus sign
                 line = line.replace("-", " -")
                 #line = line.replace("-", "-")