Energy bias fix, change the noise calculation and added a shift

pyV2DL3/vegas/fillEVENTS_not_safe.py

@@ -38,10 +38,10 @@ def __fillEVENTS_not_safe__(
     # Load header ,array info and selected event tree ( vegas > v2.5.7)
     runHeader = vegasFileIO.loadTheRunHeader()
     selectedEventsTree = vegasFileIO.loadTheCutEventTree()
+    cuts = vegasFileIO.loadTheCutsInfo()
+    arrayInfo = vegasFileIO.loadTheArrayInfo()
     qStatsData = vegasFileIO.loadTheQStatsData()
     pixelData = vegasFileIO.loadThePixelStatusData()
-    arrayInfo = vegasFileIO.loadTheArrayInfo(0)
-    cuts = vegasFileIO.loadTheCutsInfo()
 
     # Calculate time references
     startTime = runHeader.getStartTime()
@@ -84,6 +84,7 @@ def __fillEVENTS_not_safe__(
         "altArr": [],
         "energyArr": [],
         "nTelArr": [],
+        "fNoise": [],
     }
 
     # Arrays exclusive to event class mode
@@ -97,7 +98,6 @@ def __fillEVENTS_not_safe__(
         event_groups.append(deepcopy(event_arrays))
 
     logger.debug("Start filling events ...")
-
     for ev in selectedEventsTree:
         # Event reconstruction method
         if reco_type == 1:
@@ -181,6 +181,17 @@ def __fillEVENTS_not_safe__(
         this_event_group["energyArr"].append(reco.fEnergy_GeV / 1000.0)
         this_event_group["nTelArr"].append(reco.fImages)
 
+        avNoise = 0
+        nTels = 0
+        for telID in decodeConfigMask(runHeader.fRunInfo.fConfigMask):
+            avNoise += qStatsData.getCameraAverageTraceVar(
+                telID - 1, windowSizeForNoise, reco.fTime, pixelData, arrayInfo
+            )
+            nTels += 1
+
+        avNoise /= nTels
+        this_event_group["fNoise"].append(avNoise)
+
     avAlt = np.mean(avAlt)
     # Calculate average azimuth angle from average vector on a circle
     # https://en.wikipedia.org/wiki/Mean_of_circular_quantities
@@ -247,35 +258,29 @@ def __fillEVENTS_not_safe__(
         evt_dict["DEC_OBJ"] = np.rad2deg(runHeader.getSourceDec())
         evt_dict["TELLIST"] = produceTelList(runHeader.fRunInfo.fConfigMask)
         evt_dict["N_TELS"] = runHeader.pfRunDetails.fTels
-        if event_class_mode:
-            evt_dict["ENERGY"] = arr_dict["energyArr"]
 
-    avNoise = 0
-    nTels = 0
-    for telID in decodeConfigMask(runHeader.fRunInfo.fConfigMask):
-        avNoise += qStatsData.getCameraAverageTraceVarTimeIndpt(
-            telID - 1, windowSizeForNoise, pixelData, arrayInfo
-        )
-        nTels += 1
-
-    avNoise /= nTels
-    if corr_EB:
-        offset = SkyCoord(avRA * units.deg, avDec * units.deg).separation(
-            SkyCoord(arr_dict["raArr"] * units.deg, arr_dict["decArr"] * units.deg)
-        )
-        offset = offset.degree
-        evt_dict["ENERGY"] = energyBiasCorr(
-            arr_dict["energyArr"],
-            arr_dict["azArr"],
-            (90.0 - arr_dict["altArr"]),
-            avNoise,
-            offset,
-            effective_area_files[event_class_idx],
-            irf_to_store,
-            psf_king_params,
-        )
-    else:
-        evt_dict["ENERGY"] = arr_dict["energyArr"]
+        if corr_EB:
+            offset = SkyCoord(avRA * units.deg, avDec * units.deg).separation(
+                SkyCoord(arr_dict["raArr"] * units.deg, arr_dict["decArr"] * units.deg)
+            )
+            offset = offset.degree
+            #This is a problem but I don't know if it's VEGAS or V2DL3
+            #The azimuth and zenith need to be used in the previous event
+            #So this means that the first event has no reference....
+            eList=np.array([arr_dict["energyArr"][0]])
+            eList= np.append(eList,energyBiasCorr(
+                arr_dict["energyArr"],
+                arr_dict["azArr"],
+                arr_dict["altArr"],
+                arr_dict["fNoise"],
+                offset,
+                effective_area_files[event_class_idx],
+                irf_to_store,
+                psf_king_params,
+            )[1:]).flatten()
+            evt_dict["ENERGY"] = eList
+        else:
+            evt_dict["ENERGY"] = arr_dict["energyArr"]
 
     return (
         {
@@ -355,21 +360,22 @@ def energyBiasCorr(
     manager = effective_area_file.manager
     energyCorr = []
 
-    for i in range(len(energy)):
+    for i in range(0, len(energy)):
+        shift = i - 1  #this is the kludge I don't know why this needs to be done it does not make sense
         effectiveAreaParameters = ROOT.VAEASimpleParameterData()
-        effectiveAreaParameters.fAzimuth = azimuth[i]
-        effectiveAreaParameters.fZenith = zenith[i]
-        effectiveAreaParameters.fNoise = noise
-        effectiveAreaParameters.fOffset = offset[i]
+        effectiveAreaParameters.fAzimuth = azimuth[shift]
+        effectiveAreaParameters.fZenith = 90 - zenith[shift]
+        effectiveAreaParameters.fNoise = noise[shift]
+        #effectiveAreaParameters.fOffset = offset[i]
         effectiveAreaParameters = manager.getVectorParamsFromSimpleParameterData(
             effectiveAreaParameters
         )
+
         correction = manager.getCorrectionForExperimentalBias(
             effectiveAreaParameters, energy[i] * 1000
         )
+        logger.debug(f"Correction = {correction}, Original E = {energy[i]*1000}")
         energytemp = energy[i] / correction
-        logger.debug(energy[i])
-        logger.debug(str(energytemp))
         energyCorr.append(energytemp)
 
     return energyCorr