Containers for waveforms and charge extraction

.gitignore

@@ -86,4 +86,11 @@ target
 notebooks/*.html
 notebooks/*.png
 
+nectarchain/calibration/test
+**.log
+**.pdf
+**.csv
+
+#VScode
+.vscode/
 

nectarchain/calibration/README.md

@@ -0,0 +1,7 @@
+Environment variable needed: 
+export NECTARCAMDATA="path to local NectarCam data, it can contain fits.fz run files, WaveformsContainer or ChargeContainer FITS files"
+
+Environment variables which can be defined: 
+export NECTARCHAIN_TEST="path to test for nectarchain"
+export NECTARCHAIN_LOG="path to log for nectarchain"
+export NECTARCHAIN_FIGURES="path to log figures for nectarchain"

nectarchain/calibration/container/__init__.py

@@ -0,0 +1,2 @@
+from .waveforms import *
+from .charge import *

nectarchain/calibration/container/charge_extractor.py

@@ -0,0 +1,95 @@
+from ctapipe.image import ImageExtractor
+import numpy as np
+import logging
+from scipy.interpolate import InterpolatedUnivariateSpline
+from scipy.signal import find_peaks
+
+from numba import guvectorize
+
+logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s')
+log = logging.getLogger(__name__)
+log.handlers = logging.getLogger('__main__').handlers
+
+__all__ = ['gradient_extractor']
+
+class gradient_extractor(ImageExtractor) : 
+
+
+    fixed_window=np.uint8(16)
+    height_peak=np.uint8(10)
+
+    def __call__(self, waveforms, telid, selected_gain_channel,substract_ped = False):
+
+        shape = waveforms.shape
+        waveforms = waveforms.reshape(shape[0]*shape[1],shape[2])
+
+        if substract_ped : 
+            log.info('substracting pedestal')
+            #calculate pedestal
+            ped_mean = np.mean(waveforms[:,:16],axis = 1).reshape(shape[0]*shape[1],1)
+
+            y = waveforms - (ped_mean)@(np.ones(1,shape[2])) # waveforms without pedestal
+        else : 
+            log.info('do not substract pedestal')
+            y = waveforms
+
+        waveforms.reshape(shape[0],shape[1],shape[2])
+
+        peak_time, charge_integral, charge_sum, charge_integral_fixed_window, charge_sum_fixed_window = extract_charge(waveforms, self.height_peak, self.fixed_window)
+
+        return peak_time, charge_integral
+
+
+
+
+@guvectorize(
+    [
+        (np.uint16[:], np.uint8, np.uint8, np.uint16[:]),
+    ],
+    "(n,p,s),(),()->(n,p)",
+    nopython=True,
+    cache=True,
+)
+
+def extract_charge(y,height_peak,fixed_window) : 
+    x = np.linspace(0, len(y), len(y))
+    xi = np.linspace(0, len(y), 251)
+    ius = InterpolatedUnivariateSpline(x, y)
+    yi = ius(xi)
+    peaks, _ = find_peaks(yi, height=height_peak)
+    # find the peak
+    if len(peaks)>0:
+        # find the max peak
+        # max_peak = max(yi[peaks])
+        max_peak_index = np.argmax(yi[peaks], axis=0)
+        # Check if there is not a peak but a plateaux
+        # 1. divide for the maximum to round to the first digit
+        # 2. create a new array with normalized and rounded values, yi_rounded
+        yi_rounded = np.around(yi[peaks]/max(yi[peaks]),1)
+        maxima_peak_index = np.argwhere(yi_rounded == np.amax(yi_rounded))
+        if len(maxima_peak_index)>1:
+            # saturated event
+            max_peak_index = int(np.median(maxima_peak_index))
+            # width_peak = 20
+        if (xi[peaks[max_peak_index]]>20)&(xi[peaks[max_peak_index]]<40):                                                                     # Search the adaptive integration window
+            # calculate total gradients (not used, only for plot)
+            yi_grad_tot = np.gradient(yi, 1)
+            maxposition = peaks[max_peak_index]
+            # calcualte grandients starting from the max peak and going to the left to find the left margin of the window
+            yi_left = yi[:maxposition]
+            yi_grad_left = np.gradient(yi_left[::-1], 0.9)
+            change_grad_pos_left = (np.where(yi_grad_left[:-1] * yi_grad_left[1:] < 0 )[0] +1)[0]
+            # calcualte grandients starting from the max peak and going to the right to find the right margin of the window
+            yi_right = yi[maxposition:]
+            yi_grad_right = np.gradient(yi_right, 0.5)
+            change_grad_pos_right = (np.where(yi_grad_right[:-1] * yi_grad_right[1:] < 0 )[0] +1)[0]
+            charge_integral = ius.integral(xi[peaks[max_peak_index]-(change_grad_pos_left)], xi[peaks[max_peak_index]+ change_grad_pos_right])
+            charge_sum = yi[(peaks[max_peak_index]-(change_grad_pos_left)):peaks[max_peak_index]+change_grad_pos_right].sum()/(change_grad_pos_left+change_grad_pos_right)    # simple sum integration
+            adaptive_window = change_grad_pos_right + change_grad_pos_left
+            window_right = (fixed_window-6)/2
+            window_left = (fixed_window-6)/2 + 6
+            charge_integral_fixed_window = ius.integral(xi[peaks[max_peak_index]-(window_left)], xi[peaks[max_peak_index]+window_right])
+            charge_sum_fixed_window = yi[(peaks[max_peak_index]-(window_left)):peaks[max_peak_index]+window_right].sum()/(fixed_window)    # simple sum integration
+    else:
+        log.info('No peak found, maybe it is a pedestal or noisy run!')
+    return adaptive_window, charge_integral, charge_sum, charge_integral_fixed_window, charge_sum_fixed_window

nectarchain/calibration/container/utils.py

@@ -0,0 +1,177 @@
+from email.generator import Generator
+import os
+import glob
+import re
+import mechanize
+import requests
+import browser_cookie3
+
+from DIRAC.Interfaces.API.Dirac import Dirac
+from pathlib import Path
+from typing import List,Tuple
+
+
+import logging
+logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s')
+log = logging.getLogger(__name__)
+log.handlers = logging.getLogger('__main__').handlers
+
+__all__ = ['DataManagment','ChainGenerator']
+
+class DataManagment() :
+    @staticmethod
+    def findrun(run_number : int,search_on_GRID = True) -> Tuple[Path,List[Path]]: 
+        """method to find in NECTARCAMDATA the list of *.fits.fz files associated to run_number
+
+        Args:
+            run_number (int): the run number
+
+        Returns:
+            (PosixPath,list): the path list of *fits.fz files
+        """
+        basepath=os.environ['NECTARCAMDATA']
+        list = glob.glob(basepath+'**/*'+str(run_number)+'*.fits.fz',recursive=True)
+        list_path = [Path(chemin) for chemin in list]
+        if len(list_path) == 0 : 
+            e = FileNotFoundError(f"run {run_number} is not present in {basepath}")
+            if search_on_GRID : 
+                log.warning(e,exc_info=True)
+                log.info('will search files on GRID and fetch them')
+                lfns = DataManagment.get_GRID_location(run_number)
+                DataManagment.getRunFromDIRAC(lfns)
+                list = glob.glob(basepath+'**/*'+str(run_number)+'*.fits.fz',recursive=True)
+                list_path = [Path(chemin) for chemin in list]
+            else : 
+                log.error(e,exc_info=True)
+                raise e
+
+
+
+        name = list_path[0].name.split(".")
+        name[2] = "*"
+        name = Path(str(list_path[0].parent))/(f"{name[0]}.{name[1]}.{name[2]}.{name[3]}.{name[4]}")
+        log.info(f"Found {len(list_path)} files matching {name}")
+
+        #to sort list path
+        _sorted = sorted([[file,int(file.suffixes[1][1:])] for file in list_path])
+        list_path = [_sorted[i][0] for i in range(len(_sorted))]
+
+        return name,list_path
+
+    @staticmethod
+    def getRunFromDIRAC(lfns : list): 
+        """method do get run files from GRID-EGI from input lfns
+
+        Args:
+            lfns (list): list of lfns path
+        """
+
+        dirac = Dirac()
+        for lfn in lfns :
+            if not(os.path.exists(f'{os.environ["NECTARCAMDATA"]}/{os.path.basename(lfn)}')):
+                dirac.getFile(lfn=lfn,destDir=os.environ["NECTARCAMDATA"],printOutput=True)
+
+
+    @staticmethod
+    def get_GRID_location(run_number : int,output_lfns = True, username = None,password = None) : 
+        """method to get run location on GRID from Elog (work in progress!)
+
+        Args:
+            run_number (int): run number
+            output_lfns (bool, optional): if True, return lfns path of fits.gz files, else return parent directory of run location. Defaults to True.
+            username (_type_, optional): username for Elog login. Defaults to None.
+            password (_type_, optional): password for Elog login. Defaults to None.
+
+        Returns:
+            _type_: _description_
+        """
+
+        url = "http://nectarcam.in2p3.fr/elog/nectarcam-data-qm/?cmd=Find"
+
+        url_run = f"http://nectarcam.in2p3.fr/elog/nectarcam-data-qm/?mode=full&reverse=0&reverse=1&npp=20&subtext=%23{run_number}"
+
+        #try to acces data by getting cookies from firefox and Chrome
+        log.debug('try to get data with cookies from Firefox abnd Chrome')
+        cookies = browser_cookie3.load()
+        req = requests.get(f'http://nectarcam.in2p3.fr/elog/nectarcam-data-qm/?jcmd=&mode=Raw&attach=1&printable=1&reverse=0&reverse=1&npp=20&ma=&da=&ya=&ha=&na=&ca=&last=&mb=&db=&yb=&hb=&nb=&cb=&Author=&Setup=&Category=&Keyword=&Subject=&ModuleCount=&subtext=%23{run_number}',cookies = cookies)
+
+        if "<title>ELOG Login</title>" in req.text : 
+            log.debug('log to Elog with cookies impossible, try again with password')
+            #log to Elog
+            br = mechanize.Browser()
+            br.open(url)
+
+            form = br.select_form('form1')
+
+            for i in range(4) : 
+                log.debug(br.form.find_control(nr=i).name)
+
+            br.form['uname'] = username
+            br.form['upassword'] = password
+            br.method = "POST"
+
+            req = br.submit()
+            #html_page = req.get_data()
+            cookies = br._ua_handlers['_cookies'].cookiejar
+
+            #get data
+            req = requests.get(f'http://nectarcam.in2p3.fr/elog/nectarcam-data-qm/?jcmd=&mode=Raw&attach=1&printable=1&reverse=0&reverse=1&npp=20&ma=&da=&ya=&ha=&na=&ca=&last=&mb=&db=&yb=&hb=&nb=&cb=&Author=&Setup=&Category=&Keyword=&Subject=&ModuleCount=&subtext=%23{run_number}',cookies = cookies)
+
+
+        lines = req.text.split('\r\n')
+
+        url_data = None
+        for line in lines : 
+            if '<p>' in line : 
+                url_data = line.split("</p>")[0].split('FC:')[1]
+                break
+
+        if output_lfns : 
+            try : 
+                #Dirac
+                dirac = Dirac()
+                loc = f"/vo.cta.in2p3.fr/nectarcam/{url_data.split('/')[-2]}/{url_data.split('/')[-1]}"
+                res = dirac.listCatalogDirectory(loc, printOutput=True)
+
+                lfns = []
+                for key in res['Value']['Successful'][loc]['Files'].keys():
+                    if str(run_number) in key and "fits.fz" in key : 
+                        lfns.append(key)
+            except Exception as e : 
+                log.error(e,exc_info = True)
+            return lfns
+        else : 
+            return url_data
+
+
+class ChainGenerator():
+    @staticmethod
+    def chain(a : Generator ,b : Generator) :
+        """generic metghod to chain 2 generators
+
+        Args:
+            a (Generator): generator to chain
+            b (Generator): generator to chain
+
+        Yields:
+            Generator: a chain of a and b
+        """
+        yield from a
+        yield from b
+
+
+    @staticmethod
+    def chainEventSource(list : list,max_events : int = None) : #useless with ctapipe_io_nectarcam.NectarCAMEventSource
+        """recursive method to chain a list of ctapipe.io.EventSource (which may be associated to the list of *.fits.fz file of one run)
+
+        Args:
+            list (EventSource): a list of EventSource
+
+        Returns:
+            Generator: a generator which chains EventSource
+        """
+        if len(list) == 2 :
+            return ChainGenerator.chain(list[0],list[1])
+        else :
+            return ChainGenerator.chain(list[0],ChainGenerator.chainEventSource(list[1:]))
+