Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Gain computation with SPE fit and Photo-statistic methods #38

Merged
merged 100 commits into from
Mar 17, 2023
Merged
Changes from 97 commits
Commits
Show all changes
100 commits
Select commit Hold shift + click to select a range
765a802
-NectarGAIN
guillaumegrolleron Nov 9, 2022
550f7a7
-clean
guillaumegrolleron Nov 9, 2022
269f6b7
cleaning
guillaumegrolleron Dec 14, 2022
8fa7805
-new code for the SPE fit
guillaumegrolleron Dec 14, 2022
72f3cd2
-utils for NectarGain
guillaumegrolleron Dec 14, 2022
d952ce1
-implementation of photo statistic method
guillaumegrolleron Dec 14, 2022
c0336d3
-NectarGAIN
guillaumegrolleron Nov 9, 2022
aabe571
-typo log
guillaumegrolleron Dec 15, 2022
59dd96c
-bugfix gain computation
guillaumegrolleron Dec 16, 2022
ff8951e
-add parent class for SPE fit
guillaumegrolleron Dec 16, 2022
b66ff99
-update to combined SPE fit
guillaumegrolleron Dec 16, 2022
2f29b77
-add not normalized gaussian
guillaumegrolleron Jan 3, 2023
362f536
-parent class NectarGainSPE
guillaumegrolleron Jan 3, 2023
a33f69f
-cleaning child methods following parent class
guillaumegrolleron Jan 3, 2023
782d8b9
NectarGainSPE fit combined for 1000V and 1400V
guillaumegrolleron Jan 3, 2023
3d55713
minor changes
guillaumegrolleron Jan 3, 2023
e4ebedd
-luminosity can now be different
guillaumegrolleron Jan 3, 2023
8dd7886
-minor change for prefit parameters initialization
guillaumegrolleron Jan 3, 2023
cb239d3
-add class for SPE fit for run at nominal
guillaumegrolleron Jan 3, 2023
381905a
-add write method in photostat
guillaumegrolleron Jan 3, 2023
1ab2945
user script modification
guillaumegrolleron Jan 4, 2023
a6269dc
update of Photostat with container
guillaumegrolleron Jan 4, 2023
c102c0f
clean
guillaumegrolleron Jan 4, 2023
834adfa
-add linear interpolation to correlation plot
guillaumegrolleron Jan 6, 2023
e048029
-treatment of broken pixels implemented
guillaumegrolleron Jan 6, 2023
b72aa3a
-error raised if SPE fit results pixels id
guillaumegrolleron Jan 12, 2023
4d9eb1f
-NectarGAIN
guillaumegrolleron Nov 9, 2022
1f0c443
-clean
guillaumegrolleron Nov 9, 2022
8cebc5f
cleaning
guillaumegrolleron Dec 14, 2022
3556893
-new code for the SPE fit
guillaumegrolleron Dec 14, 2022
0af4ca8
-utils for NectarGain
guillaumegrolleron Dec 14, 2022
5183136
-implementation of photo statistic method
guillaumegrolleron Dec 14, 2022
15c0244
-NectarGAIN
guillaumegrolleron Nov 9, 2022
cabeaa3
-typo log
guillaumegrolleron Dec 15, 2022
d7d7ebe
-bugfix gain computation
guillaumegrolleron Dec 16, 2022
547d242
-add parent class for SPE fit
guillaumegrolleron Dec 16, 2022
924de5b
-update to combined SPE fit
guillaumegrolleron Dec 16, 2022
3f95dda
-add not normalized gaussian
guillaumegrolleron Jan 3, 2023
28af7c7
-parent class NectarGainSPE
guillaumegrolleron Jan 3, 2023
08b6300
-cleaning child methods following parent class
guillaumegrolleron Jan 3, 2023
9bf3ea2
NectarGainSPE fit combined for 1000V and 1400V
guillaumegrolleron Jan 3, 2023
62bcaea
minor changes
guillaumegrolleron Jan 3, 2023
7646d03
-luminosity can now be different
guillaumegrolleron Jan 3, 2023
86d2ce0
-minor change for prefit parameters initialization
guillaumegrolleron Jan 3, 2023
967c959
-add class for SPE fit for run at nominal
guillaumegrolleron Jan 3, 2023
c3965e6
-add write method in photostat
guillaumegrolleron Jan 3, 2023
6c644a5
update of Photostat with container
guillaumegrolleron Jan 4, 2023
a40a3a9
clean
guillaumegrolleron Jan 4, 2023
eb11174
-add linear interpolation to correlation plot
guillaumegrolleron Jan 6, 2023
999e7c2
-treatment of broken pixels implemented
guillaumegrolleron Jan 6, 2023
0f5cfaa
-error raised if SPE fit results pixels id
guillaumegrolleron Jan 12, 2023
95a962d
-add mask to broken pixels
guillaumegrolleron Jan 18, 2023
2a50b98
-add 'is_valid' column in output table
guillaumegrolleron Jan 18, 2023
91ac964
user script modif
guillaumegrolleron Jan 23, 2023
dc9bb13
Merge branch 'SPEfit' of github.com:guillaumegrolleron/nectarchain in…
Jan 26, 2023
bf9a1f0
-first implementation of multiprocessing
Jan 26, 2023
a8ece91
-multiprocessing with pool
Jan 26, 2023
125458d
-follow up implemetation multiprocessing
Jan 30, 2023
7559feb
script for 1400V data SPE fit
Feb 1, 2023
48bb97b
-add docstring
Feb 3, 2023
aaa5a52
-multiprocessing implementation for SPE fit
Feb 3, 2023
8e34a02
-save memory in wfs computation script
Feb 3, 2023
3192dda
-improve script for gain computation
Feb 3, 2023
9bd2581
-add script for gain computation of nominal run by fixing shared
Feb 3, 2023
b6b7eb6
-bugfix : hg -> lg in histo computation
Feb 3, 2023
dd0fa3f
Merge branch 'container' into SPEfit
Feb 3, 2023
5be9ae2
improve logging get DIRAC path
Feb 3, 2023
2189e3b
Merge branch 'container' into SPEfit
Feb 4, 2023
a63e176
bugfic when updating mean starting value
Feb 5, 2023
a7c4dcc
script for gain computation
Feb 5, 2023
1803269
add pandas to dependencies
Feb 5, 2023
b96a353
-quite minuit in multiprocessing
Feb 5, 2023
9cfce1c
-cleaning
Feb 6, 2023
a55b080
-renaming scripts
Feb 6, 2023
5e2a426
bugfix masked columns in output table
Feb 6, 2023
7d62219
improve logging
Feb 6, 2023
69082c9
-bugfix : output pixels_id if SPE fit has failed
Feb 7, 2023
ddf2801
-use runs of different shapes in PhotoStat
Feb 7, 2023
058f7ea
-update logging in gain photo-stat script
Feb 7, 2023
411a594
Merge branch 'master' into SPEfit
Feb 10, 2023
34e442b
-typo
Feb 10, 2023
1b24a29
-moving in src
Feb 10, 2023
1b72d67
-clean imports
Feb 10, 2023
78a31c6
typo bugfix
Feb 10, 2023
dbc3336
PR corrections
Feb 15, 2023
dc371ef
-add docstring to link with paper Caroff et al 2019
Feb 16, 2023
b14ddfe
-corrections
Feb 16, 2023
7a5ac7c
bugfix : output path was not well built from
Feb 20, 2023
6938b27
bugfix :
Feb 20, 2023
3b3119d
more generic test path for user_scripts
Feb 21, 2023
ee1923e
new pp and n fixed value for SPE fit based on full free parameters
Feb 23, 2023
18d9049
-add parameters for charge
Feb 23, 2023
89d12cc
typo correction in Caroff et al 2019
Feb 23, 2023
68aadfb
bugfix : use correction coef between FF and Ped
Feb 23, 2023
b2886ba
charge have to be compute before outside
Feb 24, 2023
3b0551a
improvment scripts
Feb 24, 2023
9cc6011
typo
Feb 24, 2023
8f9fb1c
final correction of typo in equations
Mar 16, 2023
1a85aac
likelihood value and pvalue of fit in output
Mar 16, 2023
4e12a30
Merge branch 'master' into SPEfit
Mar 17, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .gitignore
Original file line number Diff line number Diff line change
@@ -87,7 +87,7 @@ target
notebooks/*.html
notebooks/*.png

nectarchain/calibration/test
src/nectarchain/user_scripts/**/test
**.log
**.pdf
**.csv
2 changes: 1 addition & 1 deletion environment.yml
Original file line number Diff line number Diff line change
@@ -11,7 +11,7 @@ dependencies:
- pytables>=3.7
- pytest-runner
- pip
- numba
- pandas
- pip:
- git+https://github.com/cta-observatory/ctapipe_io_nectarcam.git
- mechanize
21,066 changes: 0 additions & 21,066 deletions src/nectarchain/calibration/Illuminated_pixels_charge_ped_first_post.csv

This file was deleted.

281 changes: 281 additions & 0 deletions src/nectarchain/calibration/NectarGain/PhotoStat/PhotoStat.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,281 @@
import numpy as np
from scipy.stats import linregress
from matplotlib import pyplot as plt
import astropy.units as u
from astropy.visualization import quantity_support
from astropy.table import QTable,Column
import os
from datetime import date
from pathlib import Path

import logging
logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s')
log = logging.getLogger(__name__)
log.handlers = logging.getLogger('__main__').handlers

from ctapipe_io_nectarcam import constants

from ...container import ChargeContainer,WaveformsContainer

from ..utils.error import DifferentPixelsID

from abc import ABC

__all__ = ["PhotoStatGainFFandPed"]

class PhotoStatGain(ABC):

def _readFF(self,FFRun,maxevents: int = None,**kwargs) :
log.info('reading FF data')
method = kwargs.get('method','std')
FFchargeExtractorWindowLength = kwargs.get('FFchargeExtractorWindowLength',None)
if method != 'std' :
if FFchargeExtractorWindowLength is None :
e = Exception(f"we have to specify FFchargeExtractorWindowLength argument if charge extractor method is not FullwaveformSum")
log.error(e,exc_info=True)
raise e
else :
self.__coefCharge_FF_Ped = FFchargeExtractorWindowLength / constants.N_SAMPLES
else :
self.__coefCharge_FF_Ped = 1
if isinstance(FFRun,int) :
try :
self.FFcharge = ChargeContainer.from_file(f"{os.environ['NECTARCAMDATA']}/charges/{method}",FFRun)
log.info(f'charges have ever been computed for FF run {FFRun}')
except Exception as e :
log.error("charge have not been yet computed")
raise e

#log.info(f'loading waveforms for FF run {FFRun}')
#FFwaveforms = WaveformsContainer(FFRun,maxevents)
#FFwaveforms.load_wfs()
#if method != 'std' :
# log.info(f'computing charge for FF run {FFRun} with following method : {method}')
# self.FFcharge = ChargeContainer.from_waveforms(FFwaveforms,method = method)
#else :
# log.info(f'computing charge for FF run {FFRun} with std method')
# self.FFcharge = ChargeContainer.from_waveforms(FFwaveforms)
#log.debug('writting on disk charge for further works')
#os.makedirs(f"{os.environ['NECTARCAMDATA']}/charges/{method}",exist_ok = True)
#self.FFcharge.write(f"{os.environ['NECTARCAMDATA']}/charges/{method}",overwrite = True)

elif isinstance(FFRun,ChargeContainer):
self.FFcharge = FFRun
else :
e = TypeError("FFRun must be int or ChargeContainer")
log.error(e,exc_info = True)
raise e

def _readPed(self,PedRun,maxevents: int = None,**kwargs) :
log.info('reading Ped data')
method = 'std'#kwargs.get('method','std')
if isinstance(PedRun,int) :
try :
self.Pedcharge = ChargeContainer.from_file(f"{os.environ['NECTARCAMDATA']}/charges/{method}",PedRun)
log.info(f'charges have ever been computed for Ped run {PedRun}')
except Exception as e :
log.error("charge have not been yet computed")
raise e

#log.info(f'loading waveforms for Ped run {PedRun}')
#Pedwaveforms = WaveformsContainer(PedRun,maxevents)
#Pedwaveforms.load_wfs()
#if method != 'std' :
# log.info(f'computing charge for Ped run {PedRun} with following method : {method}')
# self.Pedcharge = ChargeContainer.from_waveforms(Pedwaveforms,method = method)
#else :
# log.info(f'computing charge for Ped run {PedRun} with std method')
# self.Pedcharge = ChargeContainer.from_waveforms(Pedwaveforms)
#log.debug('writting on disk charge for further works')
#os.makedirs(f"{os.environ['NECTARCAMDATA']}/charges/{method}",exist_ok = True)
#self.Pedcharge.write(f"{os.environ['NECTARCAMDATA']}/charges/{method}",overwrite = True)

elif isinstance(PedRun,ChargeContainer):
self.Pedcharge = PedRun
else :
e = TypeError("PedRun must be int or ChargeContainer")
log.error(e,exc_info = True)
raise e

def _readSPE(self,SPEresults) :
log.info(f'reading SPE resolution from {SPEresults}')
table = QTable.read(SPEresults)
table.sort('pixel')
self.SPEResolution = table['resolution']
self.SPEGain = table['gain']
self.SPEGain_error = table['gain_error']
self._SPEvalid = table['is_valid']
self._SPE_pixels_id = np.array(table['pixel'].value,dtype = np.uint16)

def _reshape_all(self) :
log.info("reshape of SPE, Ped and FF data with intersection of pixel ids")
FFped_intersection = np.intersect1d(self.Pedcharge.pixels_id,self.FFcharge.pixels_id)
SPEFFPed_intersection = np.intersect1d(FFped_intersection,self._SPE_pixels_id[self._SPEvalid])
self._pixels_id = SPEFFPed_intersection
log.info(f"data have {len(self._pixels_id)} pixels in common")

self._FFcharge_hg = self.FFcharge.select_charge_hg(SPEFFPed_intersection)
self._FFcharge_lg = self.FFcharge.select_charge_lg(SPEFFPed_intersection)

self._Pedcharge_hg = self.Pedcharge.select_charge_hg(SPEFFPed_intersection)
self._Pedcharge_lg = self.Pedcharge.select_charge_lg(SPEFFPed_intersection)

#self._mask_FF = np.array([self.FFcharge.pixels_id[i] in SPEFFPed_intersection for i in range(self.FFcharge.npixels)],dtype = bool)
#self._mask_Ped = np.array([self.Pedcharge.pixels_id[i] in SPEFFPed_intersection for i in range(self.Pedcharge.npixels)],dtype = bool)
self._mask_SPE = np.array([self._SPE_pixels_id[i] in SPEFFPed_intersection for i in range(len(self._SPE_pixels_id))],dtype = bool)



def create_output_table(self) :
self._output_table = QTable()
self._output_table.meta['npixel'] = self.npixels
self._output_table.meta['comments'] = f'Produced with NectarGain, Credit : CTA NectarCam {date.today().strftime("%B %d, %Y")}'

self._output_table.add_column(Column(np.ones((self.npixels),dtype = bool),"is_valid",unit = u.dimensionless_unscaled))
self._output_table.add_column(Column(self.pixels_id,"pixel",unit = u.dimensionless_unscaled))
self._output_table.add_column(Column(np.empty((self.npixels),dtype = np.float64),"high gain",unit = u.dimensionless_unscaled))
self._output_table.add_column(Column(np.empty((self.npixels,2),dtype = np.float64),"high gain error",unit = u.dimensionless_unscaled))
self._output_table.add_column(Column(np.empty((self.npixels),dtype = np.float64),"low gain",unit = u.dimensionless_unscaled))
self._output_table.add_column(Column(np.empty((self.npixels,2),dtype = np.float64),"low gain error",unit = u.dimensionless_unscaled))

def run(self,**kwargs):
log.info('running photo statistic method')

self._output_table["high gain"] = self.gainHG
self._output_table["low gain"] = self.gainLG
#self._output_table["is_valid"] = self._SPEvalid

def save(self,path,**kwargs) :
path = Path(path)
os.makedirs(path,exist_ok = True)
log.info(f'data saved in {path}')
self._output_table.write(f"{path}/output_table.ecsv", format='ascii.ecsv',overwrite = kwargs.get("overwrite",False))

def plot_correlation(self) :
mask = (self._output_table["high gain"]>20) * (self.SPEGain[self._mask_SPE]>0) * (self._output_table["high gain"]<80) * self._output_table['is_valid']
a, b, r, p_value, std_err = linregress(self._output_table["high gain"][mask], self.SPEGain[self._mask_SPE][mask],'greater')
x = np.linspace(self._output_table["high gain"][mask].min(),self._output_table["high gain"][mask].max(),1000)
y = lambda x: a * x + b
with quantity_support() :
fig,ax = plt.subplots(1,1,figsize=(8, 6))
ax.scatter(self._output_table["high gain"][mask],self.SPEGain[self._mask_SPE][mask],marker =".")
ax.plot(x,y(x),color = 'red', label = f"linear fit,\n a = {a:.2e},\n b = {b:.2e},\n r = {r:.2e},\n p_value = {p_value:.2e},\n std_err = {std_err:.2e}")
ax.plot(x,x,color = 'black',label = "y = x")
ax.set_xlabel("Gain Photo stat (ADC)", size=15)
ax.set_ylabel("Gain SPE fit (ADC)", size=15)
#ax.set_xlim(xmin = 0)
#ax.set_ylim(ymin = 0)

ax.legend(fontsize=15)
return fig

@property
def npixels(self) : return self._pixels_id.shape[0]

@property
def pixels_id(self) : return self._pixels_id

@property
def sigmaPedHG(self) : return np.std(self._Pedcharge_hg ,axis = 0) * np.sqrt(self.__coefCharge_FF_Ped)

@property
def sigmaChargeHG(self) : return np.std(self._FFcharge_hg - self.meanPedHG, axis = 0)

@property
def meanPedHG(self) : return np.mean(self._Pedcharge_hg ,axis = 0) * self.__coefCharge_FF_Ped

@property
def meanChargeHG(self) : return np.mean(self._FFcharge_hg - self.meanPedHG, axis = 0)

@property
def BHG(self) :
min_events = np.min((self._FFcharge_hg.shape[0],self._Pedcharge_hg.shape[0]))
upper = (np.power(self._FFcharge_hg.mean(axis = 1)[:min_events] - self._Pedcharge_hg.mean(axis = 1)[:min_events] * self.__coefCharge_FF_Ped - self.meanChargeHG.mean(),2)).mean(axis = 0)
lower = np.power(self.meanChargeHG.mean(),2)#np.power(self.meanChargeHG,2)#np.power(self.meanChargeHG.mean(),2)
return np.sqrt(upper/lower)

@property
def gainHG(self) :
return ((np.power(self.sigmaChargeHG,2) - np.power(self.sigmaPedHG,2) - np.power(self.BHG * self.meanChargeHG,2))
/(self.meanChargeHG * (1 + np.power(self.SPEResolution[self._mask_SPE],2))))


@property
def sigmaPedLG(self) : return np.std(self._Pedcharge_lg ,axis = 0) * np.sqrt(self.__coefCharge_FF_Ped)

@property
def sigmaChargeLG(self) : return np.std(self._FFcharge_lg - self.meanPedLG,axis = 0)

@property
def meanPedLG(self) : return np.mean(self._Pedcharge_lg,axis = 0) * self.__coefCharge_FF_Ped

@property
def meanChargeLG(self) : return np.mean(self._FFcharge_lg - self.meanPedLG,axis = 0)

@property
def BLG(self) :
min_events = np.min((self._FFcharge_lg.shape[0],self._Pedcharge_lg.shape[0]))
upper = (np.power(self._FFcharge_lg.mean(axis = 1)[:min_events] - self._Pedcharge_lg.mean(axis = 1)[:min_events] * self.__coefCharge_FF_Ped - self.meanChargeLG.mean(),2)).mean(axis = 0)
lower = np.power(self.meanChargeLG.mean(),2) #np.power(self.meanChargeLG,2) #np.power(self.meanChargeLG.mean(),2)
return np.sqrt(upper/lower)

@property
def gainLG(self) : return ((np.power(self.sigmaChargeLG,2) - np.power(self.sigmaPedLG,2) - np.power(self.BLG * self.meanChargeLG,2))
/(self.meanChargeLG * (1 + np.power(self.SPEResolution[self._mask_SPE],2))))



class PhotoStatGainFFandPed(PhotoStatGain):
def __init__(self, FFRun, PedRun, SPEresults : str, maxevents : int = None, **kwargs) :
self._readFF(FFRun,maxevents,**kwargs)
self._readPed(PedRun,maxevents,**kwargs)

"""
if self.FFcharge.charge_hg.shape[1] != self.Pedcharge.charge_hg.shape[1] :
e = Exception("Ped run and FF run must have the same number of pixels")
log.error(e,exc_info = True)
raise e
"""

self._readSPE(SPEresults)
##need to implement reshape of SPE results with FF and Ped pixels ids
self._reshape_all()

"""
if (self.FFcharge.pixels_id.shape[0] != self._SPE_pixels_id.shape[0]) :
e = Exception("Ped run and FF run must have the same number of pixels as SPE fit results")
log.error(e,exc_info = True)
raise e

if (self.FFcharge.pixels_id != self.Pedcharge.pixels_id).any() or (self.FFcharge.pixels_id != self._SPE_pixels_id).any() or (self.Pedcharge.pixels_id != self._SPE_pixels_id).any() :
e = DifferentPixelsID("Ped run, FF run and SPE run need to have same pixels id")
log.error(e,exc_info = True)
raise e
else :
self._pixels_id = self.FFcharge.pixels_id
"""
self.create_output_table()



class PhotoStatGainFF(PhotoStatGain):
def __init__(self, FFRun, PedRun, SPEresults : str, maxevents : int = None, **kwargs) :
e = NotImplementedError("PhotoStatGainFF is not yet implemented")
log.error(e, exc_info = True)
raise e
self._readFF(FFRun,maxevents,**kwargs)

self._readSPE(SPEresults)

"""
if self.FFcharge.pixels_id != self._SPE_pixels_id :
e = DifferentPixelsID("Ped run, FF run and SPE run need to have same pixels id")
log.error(e,exc_info = True)
raise e
else :
self._pixels_id = self.FFcharge.pixels_id
"""

self._reshape_all()

self.create_output_table()
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
from .PhotoStat import *
Original file line number Diff line number Diff line change
@@ -7,33 +7,6 @@
from iminuit import Minuit
guillaumegrolleron marked this conversation as resolved.
Show resolved Hide resolved
import random

def make_minuit_par_kwargs(parameters,parnames,LimitLow,LimitUp):
"""Create *Parameter Keyword Arguments* for the `Minuit` constructor.

updated for Minuit >2.0
"""
names = parnames
kwargs = {"names": names,"values" : {}}

for i in range(len(parameters)):
kwargs["values"][parnames[i]] = parameters[i]
min_ = None if np.isnan(LimitLow[i]) else LimitLow[i]
max_ = None if np.isnan(LimitUp[i]) else LimitUp[i]
kwargs[f"limit_{names[i]}"] = (min_, max_)
kwargs[f"error_{names[i]}"] = 0.1

return kwargs

def set_minuit_parameters_limits_and_errors(m : Minuit,parameters : dict) :
"""function to set minuit parameter limits and errors with Minuit >2.0

Args:
m (Minuit): a Minuit instance
parameters (dict): dict containing parameters names, limits errors and values
"""
for name in parameters["names"] :
m.limits[name] = parameters[f"limit_{name}"]
m.errors[name] = parameters[f"error_{name}"]



Loading