clean commit history :

-update SPE fit parameters initialization following updates of the ctapipe extractors -improve the limits of generated charge histogram plots
cta-observatory · Jul 12, 2023 · 5c05525 · 5c05525
1 parent 9b4fc9e
commit 5c05525
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 17 deletions.
diff --git a/src/nectarchain/calibration/NectarGain/SPEfit/NectarGainSPE.py b/src/nectarchain/calibration/NectarGain/SPEfit/NectarGainSPE.py
@@ -30,6 +30,9 @@
 
 class NectarGainSPE(ABC) :
     _Ncall = 4000000
+    _Windows_lenght = 40
+    _Order = 2
+
     def __init__(self) :
         #set parameters value for fit
         self.__parameters = Parameters()
@@ -104,8 +107,8 @@ def _get_mean_gaussian_fit(charge_in, histo_in ,extension = ""):
         charge = charge_in.data[~histo_in.mask]
         histo = histo_in.data[~histo_in.mask]
 
-        windows_lenght = 80
-        order = 2
+        windows_lenght = NectarGainSPE._Windows_lenght
+        order = NectarGainSPE._Order
         histo_smoothed = savgol_filter(histo, windows_lenght, order)
 
         peaks = find_peaks(histo_smoothed,10)
@@ -126,15 +129,16 @@ def _get_mean_gaussian_fit(charge_in, histo_in ,extension = ""):
 
         if log.getEffectiveLevel() == logging.DEBUG :
             log.debug('plotting figures with prefit parameters computation') 
-            fig,ax = plt.subplots(1,1,figsize = (8,8))
+            fig,ax = plt.subplots(1,1,figsize = (5,5))
             ax.errorbar(charge,histo,np.sqrt(histo),zorder=0,fmt=".",label = "data")
             ax.plot(charge,histo_smoothed,label = f'smoothed data with savgol filter (windows lenght : {windows_lenght}, order : {order})')
             ax.plot(charge,weight_gaussian(charge,coeff_mean[0],coeff_mean[1],coeff_mean[2]),label = 'gaussian fit of the SPE')
             ax.vlines(coeff_mean[1],0,peak_value,label = f'mean initial value = {coeff_mean[1] - coeff[1]:.0f}',color = "red")
             ax.add_patch(Rectangle((coeff_mean[1]-coeff_mean[2], 0), 2 * coeff_mean[2], peak_value_mean,fc=to_rgba('red', 0.5)))
+            ax.set_xlim([peak_pos - 500,None])
             ax.set_xlabel("Charge (ADC)", size=15)
             ax.set_ylabel("Events", size=15)
-            ax.legend(fontsize=15)
+            ax.legend(fontsize=7)
             os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/",exist_ok=True)
             fig.savefig(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/initialization_mean_pixel{extension}_{os.getpid()}.pdf")
             fig.clf()
@@ -149,8 +153,8 @@ def _get_pedestal_gaussian_fit(charge_in, histo_in ,extension = "") :
         charge = charge_in.data[~histo_in.mask]
         histo = histo_in.data[~histo_in.mask]
 
-        windows_lenght = 80
-        order = 2
+        windows_lenght = NectarGainSPE._Windows_lenght
+        order = NectarGainSPE._Order
         histo_smoothed = savgol_filter(histo, windows_lenght, order)
 
         peaks = find_peaks(histo_smoothed,10)
@@ -161,15 +165,16 @@ def _get_pedestal_gaussian_fit(charge_in, histo_in ,extension = "") :
 
         if log.getEffectiveLevel() == logging.DEBUG :
             log.debug('plotting figures with prefit parameters computation') 
-            fig,ax = plt.subplots(1,1,figsize = (8,8))
+            fig,ax = plt.subplots(1,1,figsize = (5,5))
             ax.errorbar(charge,histo,np.sqrt(histo),zorder=0,fmt=".",label = "data")
             ax.plot(charge,histo_smoothed,label = f'smoothed data with savgol filter (windows lenght : {windows_lenght}, order : {order})')
             ax.plot(charge,weight_gaussian(charge,coeff[0],coeff[1],coeff[2]),label = 'gaussian fit of the pedestal, left tail only')
+            ax.set_xlim([peak_pos - 500,None])
             ax.vlines(coeff[1],0,peak_value,label = f'pedestal initial value = {coeff[1]:.0f}',color = 'red')
             ax.add_patch(Rectangle((coeff[1]-coeff[2], 0), 2 * coeff[2], peak_value,fc=to_rgba('red', 0.5)))
             ax.set_xlabel("Charge (ADC)", size=15)
             ax.set_ylabel("Events", size=15)
-            ax.legend(fontsize=15)
+            ax.legend(fontsize=7)
             os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/",exist_ok=True)
             fig.savefig(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/initialization_pedestal_pixel{extension}_{os.getpid()}.pdf")
             fig.clf()

diff --git a/src/nectarchain/calibration/NectarGain/SPEfit/NectarGainSPE_singlerun.py b/src/nectarchain/calibration/NectarGain/SPEfit/NectarGainSPE_singlerun.py
@@ -4,7 +4,7 @@
 from iminuit import Minuit
 import random
 import astropy.units as u
-from astropy.table import QTable,Column,MaskedColumn
+from astropy.table import Table,QTable,Column,MaskedColumn
 import yaml
 import os
 from datetime import date
@@ -457,17 +457,18 @@ def _run_obs_static(cls,it : int, funct,parameters : Parameters, pixels_id : int
 
 
             if kwargs.get('figpath',0) != 0 :
-                fig,ax = plt.subplots(1,1,figsize=(8, 6))
+                fig,ax = plt.subplots(1,1,figsize=(8, 8))
                 ax.errorbar(charge,histo,np.sqrt(histo),zorder=0,fmt=".",label = "data")
                 ax.plot(charge,
                     np.trapz(histo,charge)*MPE2(charge,parameters['pp'].value, parameters['resolution'].value, parameters['mean'].value, parameters['n'].value, parameters['pedestal'].value, parameters['pedestalWidth'].value, parameters['luminosity'].value),
                     zorder=1,
                     linewidth=2,
-                    label = f"SPE model fit \n gain : {gain[0] - gain[1]:.2f} < {gain[0]:.2f} < {gain[0] + gain[2]:.2f} ADC/pe, pvalue = {Statistics.chi2_pvalue(ndof,fit.fval)},\n likelihood = {fit.fval:.2f}")
+                    label = f"SPE model fit \n gain : {gain[0] - gain[1]:.2f} < {gain[0]:.2f} < {gain[0] + gain[2]:.2f} ADC/pe,\n pvalue = {Statistics.chi2_pvalue(ndof,fit.fval)},\n likelihood = {fit.fval:.2f}")
                 ax.set_xlabel("Charge (ADC)", size=15)
                 ax.set_ylabel("Events", size=15)
                 ax.set_title(f"SPE fit pixel {it} with pixel_id : {pixels_id}")
-                ax.legend(fontsize=15)
+                ax.set_xlim([parameters['pedestal'].value - 6 * parameters['pedestalWidth'].value, None])
+                ax.legend(fontsize=18)
                 os.makedirs(kwargs.get('figpath'),exist_ok = True)
                 fig.savefig(f"{kwargs.get('figpath')}/fit_SPE_pixel{pixels_id}.pdf")
                 fig.clf()
@@ -503,17 +504,18 @@ def _run_obs(self,pixel,prescan = False,**kwargs) :
 
 
             if kwargs.get('figpath',0) != 0 :
-                fig,ax = plt.subplots(1,1,figsize=(8, 6))
+                fig,ax = plt.subplots(1,1,figsize=(8, 8))
                 ax.errorbar(self.charge[pixel],self.histo[pixel],np.sqrt(self.histo[pixel]),zorder=0,fmt=".",label = "data")
                 ax.plot(self.charge[pixel],
                     np.trapz(self.histo[pixel],self.charge[pixel])*MPE2(self.charge[pixel],self.__pp.value,self.__resolution.value,self.__mean.value,self.__n.value,self.pedestal.value,self.__pedestalWidth.value,self.__luminosity.value),
                     zorder=1,
                     linewidth=2,
-                    label = f"SPE model fit \n gain : {self.__gain[pixel,0] - self.__gain[pixel,1]:.2f} < {self.__gain[pixel,0]:.2f} < {self.__gain[pixel,0] + self.__gain[pixel,2]:.2f} ADC/pe, pvalue = {Statistics.chi2_pvalue(ndof,fit.fval)},\n likelihood = {fit.fval:.2f}")
+                    label = f"SPE model fit \n gain : {self.__gain[pixel,0] - self.__gain[pixel,1]:.2f} < {self.__gain[pixel,0]:.2f} < {self.__gain[pixel,0] + self.__gain[pixel,2]:.2f} ADC/pe,\n pvalue = {Statistics.chi2_pvalue(ndof,fit.fval)},\n likelihood = {fit.fval:.2f}")
                 ax.set_xlabel("Charge (ADC)", size=15)
                 ax.set_ylabel("Events", size=15)
                 ax.set_title(f"SPE fit pixel : {pixel} (pixel id : {self.pixels_id[pixel]})")
-                ax.legend(fontsize=15)
+                ax.set_xlim([self.pedestal.value - 6 * self.pedestalWidth.value, None])
+                ax.legend(fontsize=18)
                 os.makedirs(kwargs.get('figpath'),exist_ok = True)
                 fig.savefig(f"{kwargs.get('figpath')}/fit_SPE_pixel{self.pixels_id[pixel]}.pdf")
                 fig.clf()
@@ -571,7 +573,7 @@ def _update_parameters_prefit_static(cls,it : int, parameters : Parameters, char
         log.debug(f"pedestalWidth updated : {pedestalWidth}")
         try : 
             coeff,var_matrix =  NectarGainSPE._get_mean_gaussian_fit(charge,histo,f'{it}_nominal')
-            if coeff[1] - pedestal.value < 0 : raise Exception("mean gaussian fit not good")
+            if (coeff[1] - pedestal.value < 0) or ((coeff[1] - coeff[2]) - pedestal.max < 0) : raise Exception("mean gaussian fit not good")
             mean = parameters['mean']
             mean.value = coeff[1] - pedestal.value
             mean.min = (coeff[1] - coeff[2]) - pedestal.max
@@ -785,7 +787,7 @@ def run(self,pixel : int = None,multiproc = False, **kwargs):
         super().run(pixel,multiproc,**kwargs)
 
     def _run_obs(self,pixel,prescan = False,**kwargs) : 
-        if self.__nectarGainSPEresult[pixel]['is_valid'].value : 
+        if self.__nectarGainSPEresult[pixel]['is_valid'] : 
             kwargs['pixel_id'] = self.pixels_id[pixel]
             super()._run_obs(pixel,prescan,**kwargs)
         else :