Merge branch 'dev_sim2root' of https://github.com/grand-mother/grand …

…into dev_sim2root

grand/geo/coordinates.py

@@ -571,7 +571,9 @@ class Geodetic(GeodeticRepresentation):
     Longitude:	Angle east and west of the Prime Meridian. The Prime Meridian
                             is a north-south line that passes through Greenwich, UK.
                             +ve to the east of the Prime Meridian, -ve to the west.
-                            Range: -180 deg to +180 deg.
+                            Range: 0 deg to 360 deg positive or negative. 
+                            Note that coordinate transformation is possible for +ve 0 to 360 deg.
+                            So negative values are changed to positive by adding 360.
     Height:	Also called altitude or elevation, this represents the height above
                     the Earth ellipsoid, measured in meters. The Earth ellipsoid is a
                     mathematical surface defined by a semi-major axis and a semi-minor axis.
@@ -702,6 +704,10 @@ def __init__(
 
         if isinstance(latitude, (Number, np.ndarray)):
             # use setter to replace placeholder coordinates values with the real values.
+            # RK: +ve 0 to 360 were only accepted for coordinate transformation. 
+            #     Now both +ve and -ve values are accepted for longitudes.
+            longitude = 360+longitude if longitude<0 else longitude 
+
             self.latitude = latitude
             self.longitude = longitude
             self.height = height

grand/sim/detector/process_ant.py

@@ -276,6 +276,8 @@ def compute_voltage(
         # frame is shower frame. self.frame is antenna frame.
         logger.debug(f"pos {self.pos}")
         if (not np.all(np.isfinite(self.pos))) or (not np.all(np.isfinite(frame))): # which one
+            print("pos",self.pos)
+            print("frame",frame)
             raise MissingFrameError("missing antenna or shower frame")
 
         # Compute the voltage. input fft_leff and field are in shower frame.

grand/sim/efield2voltage.py

@@ -345,6 +345,8 @@ def compute_voc_du(self, du_idx):
 
         # ----- antenna responses -----
         # compute_voltage() --> return Voltage(t=t, V=volt_t)
+        print("here frame",self.evt_shower.frame)
+
         self.voc[du_idx, 0] = self.ant_leff_sn.compute_voltage(
             self.evt_shower.maximum, efield_idx, self.evt_shower.frame
         ).V

sim2root/CoREASRawRoot/CoreasToRawROOT.py

@@ -3,6 +3,7 @@
 ## by Jelena Köhler, @jelenakhlr
 
 import sys
+import numpy as np
 import os
 import glob
 import time #to get the unix timestamp
@@ -237,6 +238,15 @@ def CoreasToRawRoot(file, simID=None):
   ed_sum = energy_dep[:,9]
 
 
+  Egamma = np.sum(ed_gamma)
+  Eem_ion = np.sum(ed_em_ioniz)
+  Eem_cut = np.sum(ed_em_cut)
+
+  # calculate electromagnetic shower energy and convert to eV
+  Eem = (Egamma + Eem_ion + Eem_cut) * 1e9
+  print("Electromagnetic shower energy:", Eem)
+
+
   ##############################################
 
   EnergyInNeutrinos = 1. # placeholder
@@ -261,6 +271,9 @@ def CoreasToRawRoot(file, simID=None):
   site = read_site(inp_input)
   latitude, longitude, altitude = read_lat_long_alt(site)
 
+  # set altitude to simulations obslevel
+  altitude = CorePosition[2]
+
   ############################################################################################################################
   # Part B.I.ii: Create and fill the RAW Shower Tree
   ############################################################################################################################
@@ -282,6 +295,7 @@ def CoreasToRawRoot(file, simID=None):
   RawShower.rnd_seed = RandomSeed
 
   RawShower.energy_in_neutrinos = EnergyInNeutrinos
+  RawShower.energy_em = [Eem]
   RawShower.energy_primary = [Energy]
   RawShower.azimuth = azimuth
   RawShower.zenith = zenith

sim2root/Common/raw_root_trees.py

@@ -67,6 +67,9 @@ class RawShowerTree(MotherEventTree):
 
     ###X Primary energy (GeV) 
     energy_primary: StdVectorListDesc = field(default=StdVectorListDesc("float"))
+
+    ###X Electromagnetic energy (GeV) 
+    energy_em: StdVectorListDesc = field(default=StdVectorListDesc("float"))
 
     ### Shower azimuth (deg, CR convention)
     _azimuth: np.ndarray = field(default_factory=lambda: np.zeros(1, np.float32))

sim2root/Common/sim2root.py

@@ -11,7 +11,7 @@
 import raw_root_trees as RawTrees # this is here in Common
 import grand.manage_log as mlg
 import matplotlib.pyplot as plt
-from scipy.ndimage.interpolation import shift  #to shift the time trance for the trigger simulation
+from scipy.ndimage import shift  #to shift the time trance for the trigger simulation
 
 # specific logger definition for script because __mane__ is "__main__" !
 logger = mlg.get_logger_for_script(__file__)
@@ -560,6 +560,8 @@ def rawshower2grandroot(trawshower, gt):
     # ToDo: it should be a scalar on sim side
     gt.tshower.energy_primary = trawshower.energy_primary[0]
 
+    gt.tshower.energy_em = trawshower.energy_em[0]
+
     ### Shower azimuth (deg, CR convention)
     gt.tshower.azimuth = trawshower.azimuth