diff --git a/src/nectarchain/user_scripts/vmarandon/CalibrationData.py b/src/nectarchain/user_scripts/vmarandon/CalibrationData.py index a146ee5b..f574f8cf 100644 --- a/src/nectarchain/user_scripts/vmarandon/CalibrationData.py +++ b/src/nectarchain/user_scripts/vmarandon/CalibrationData.py @@ -13,12 +13,14 @@ class CalibrationCameraDisplay(CameraDisplay): def __init__(self,*args, **kwargs): super().__init__(*args,**kwargs) + self.clickfunc = None def set_function(self,func_name): self.clickfunc = func_name def on_pixel_clicked(self, pix_id): - self.clickfunc(pix_id) + if self.clickfunc is not None: + self.clickfunc(pix_id) class CalibInfo: @@ -156,8 +158,25 @@ def ShowPedestal(self): plt.show() - - +class XYTableDataElement(TimedInfo): + '''Class to store waveforms for each position of the XY tables''' + def __init__(self,startTime=None, endTime=None,bloc=None): + super().__init__(startTime,endTime) + self.waveforms = None + self.masks = None + self.averaged_waveform = None + self.bloc_number = bloc + +class XYTableDataIntegratedElement(TimedInfo): + '''Class to store integrated waveforms for each position of the XY tables''' + def __init__(self,startTime=None,endTime=None,bloc=None): + super().__init__(startTime,endTime) + self.pixels = None + self.times = None + self.pedestals = None + self.pedwidths = None + self.integrated = None + self.bloc_number = bloc class FlatFieldInfo(TimedInfo): diff --git a/src/nectarchain/user_scripts/vmarandon/DBHandler.py b/src/nectarchain/user_scripts/vmarandon/DBHandler.py index a468d03d..523daddb 100644 --- a/src/nectarchain/user_scripts/vmarandon/DBHandler.py +++ b/src/nectarchain/user_scripts/vmarandon/DBHandler.py @@ -96,14 +96,14 @@ def __get_selection_condition__(self,table): #df1b = pd.read_sql(f"SELECT * FROM 'monitoring_drawer_temperatures' WHERE time>datetime({int(test_time.timestamp())}, 'unixepoch')", con=sqlite3.connect(db_url)) time_cond = "" if self.time_start is not None: - if isinstance( self.time_start, datetime ): + if isinstance( self.time_start, datetime.datetime ): time_cond = f" WHERE time >= datetime({self.time_start.timestamp()}, 'unixepoch') " else: print(f"WARNING> {self.time_start} of type {type(self.time_start)} is of a non handled type ==> Won't be used (please correct code)") if self.time_end is not None: - if isinstance( self.time_end, datetime ): + if isinstance( self.time_end, datetime.datetime ): link_word = " WHERE " if not time_cond else " AND " time_cond = f" {link_word} time <= datetime({self.time_end.timestamp()}, 'unixepoch') " else: diff --git a/src/nectarchain/user_scripts/vmarandon/DataUtils.py b/src/nectarchain/user_scripts/vmarandon/DataUtils.py index 3e396956..a91a0ec1 100644 --- a/src/nectarchain/user_scripts/vmarandon/DataUtils.py +++ b/src/nectarchain/user_scripts/vmarandon/DataUtils.py @@ -108,7 +108,8 @@ def GetLongRunTimeEdges(run,path=None,event_type=None,delta_t_second=10.): #print(nEvents) data = DataReader(run,path=path) - data.Connect("trigger") + if not data.Connect("trigger"): + data = GetNectarCamEvents(run=run,path=path,applycalib=False) times_edges = list() @@ -119,26 +120,27 @@ def GetLongRunTimeEdges(run,path=None,event_type=None,delta_t_second=10.): # if there is a time gap of more than delta_t seconds, then we consider that this is the end of a data block delta_t = TimeDelta(delta_t_second,format="sec") - - for evt in tqdm(data,total=nEvents): - current_time = data.trigger.time - - if time_start is None: - time_start = current_time - - if evt.trigger.event_type == EventType.SKY_PEDESTAL: - - if previous_time is None: - previous_time = current_time + try: + for evt in tqdm(data,total=nEvents): + current_time = evt.trigger.time - if current_time - previous_time > delta_t: - #print(f"time: {time} previous time: {previous_time} delta: {(time - previous_time).to_value('s')}") - #if (previous_time - time) > delta_t: - times_edges.append( (time_start,previous_time) ) + if time_start is None: time_start = current_time + + if evt.trigger.event_type == event_type: - previous_time = current_time - + if previous_time is None: + previous_time = current_time + + if current_time - previous_time > delta_t: + #print(f"time: {time} previous time: {previous_time} delta: {(time - previous_time).to_value('s')}") + #if (previous_time - time) > delta_t: + times_edges.append( (time_start,previous_time) ) + time_start = current_time + + previous_time = current_time + except Exception as err: + print(f"Error while reading file: [{err}]") times_edges.append( (time_start,current_time) ) # write the last time #print(f"There is : {len(times_edges)} intervals") diff --git a/src/nectarchain/user_scripts/vmarandon/ExtractInformation2.py b/src/nectarchain/user_scripts/vmarandon/ExtractInformation2.py index 7cf3b18d..f00f6917 100644 --- a/src/nectarchain/user_scripts/vmarandon/ExtractInformation2.py +++ b/src/nectarchain/user_scripts/vmarandon/ExtractInformation2.py @@ -16,6 +16,7 @@ from ctapipe.containers import EventType from multiprocessing.dummy import Pool as ThreadPool + #from multiprocessing import Pool as ThreadPool @@ -33,7 +34,7 @@ # data_path = FindDataPath(run,args.dataPath) -def ExtractInformationSingleRun(run,data_path,dest_path,data_block,applycalib=True,keepR1=True,nnint=False): +def ExtractInformationSingleRun(run,data_path,dest_path,data_block,applycalib=True,keepR1=True,nnint=False,onlytrigger=False): sleep_time = random.uniform(0.,1.) #print(sleep_time) @@ -68,6 +69,9 @@ def ExtractInformationSingleRun(run,data_path,dest_path,data_block,applycalib=Tr doIntegration = nnint + if onlytrigger: + doIntegration = False + for evt in tqdm(events): #if data_block != 42 and data_block!=8: # break @@ -113,7 +117,10 @@ def ExtractInformationSingleRun(run,data_path,dest_path,data_block,applycalib=Tr for k in evt.keys(): if k not in exclusion: - dd[k].dump( copy.deepcopy(getattr(evt,k)), time=event_time ) + if onlytrigger and k!="trigger": + continue + else: + dd[k].dump( copy.deepcopy(getattr(evt,k)), time=event_time ) def TrueOrFalse(arg): ua = str(arg).upper() @@ -140,6 +147,7 @@ def ExtractInformation(arglist): p.add_argument("--keep-r1",dest='keepR1',type=str,default="True",help="Save the R1 data if True") p.add_argument("--split",dest='split',action='store_true',help='Split the files per groups. 0-1 together in a file, 2-3 in another, etc... Need the ctapipe_io_nectarcam version compatible with ctapipe 0.18') p.add_argument("--nnint",dest='nnint',action='store_true',help='Do an integration of the data using Next Neighbor Peak Search. At the moment hard coded to be 10 ns -4 and +6 ns after the max. Will create charge and TO data set') + p.add_argument("--only-trigger",dest='onlytrig',action='store_true',help='Extract only the trigger information to a file. Useful for big runs') args = p.parse_args(arglist) @@ -164,7 +172,7 @@ def ExtractInformation(arglist): dest_path = args.destPath - if args.split: + if False and args.split: runs = list() paths = list() @@ -173,6 +181,7 @@ def ExtractInformation(arglist): calib = list() keepR1 = list() nnints = list() + trigonly = list() for block in range(GetNumberOfDataBlocks(run,data_path)): #for block in range(8): @@ -184,13 +193,14 @@ def ExtractInformation(arglist): calib.append(applyCalib) keepR1.append( keepR1) nnints.append(args.nnint) + trigonly.append(args.onlytrig) # Make the Pool of workers - pool = ThreadPool(1) + pool = ThreadPool(4) # Open the URLs in their own threads # and return the results - results = pool.starmap(ExtractInformationSingleRun, zip(runs,paths,dest_paths,blocks,calib,keepR1,nnints) ) + results = pool.starmap(ExtractInformationSingleRun, zip(runs,paths,dest_paths,blocks,calib,keepR1,nnints,trigonly) ) # Close the pool and wait for the work to finish pool.close() @@ -198,7 +208,14 @@ def ExtractInformation(arglist): #ExtractInformationSingleRun(run,args.dataPath,args.destPath,data_block=block) else: - ExtractInformationSingleRun(run=run,data_path=data_path,dest_path=dest_path,data_block=-1,applycalib=applyCalib,keepR1=keepR1,nnint=args.nnint) + if args.split: + nBlocks = GetNumberOfDataBlocks(run,data_path) + for block in range(nBlocks): + print(f'block: {block+1}/{nBlocks}') + ExtractInformationSingleRun(run=run,data_path=data_path,dest_path=dest_path,data_block=block,applycalib=applyCalib,keepR1=keepR1,nnint=args.nnint,onlytrigger=args.onlytrig) + else: + ExtractInformationSingleRun(run=run,data_path=data_path,dest_path=dest_path,data_block=-1,applycalib=applyCalib,keepR1=keepR1,nnint=args.nnint,onlytrigger=args.onlytrig) + #def ExtractInformationSingleRun(run,data_path,dest_path,data_block,applycalib=True,keepR1=True): diff --git a/src/nectarchain/user_scripts/vmarandon/FitSPEFFMultiProcess.py b/src/nectarchain/user_scripts/vmarandon/FitSPEFFMultiProcess.py new file mode 100644 index 00000000..a52f3ad3 --- /dev/null +++ b/src/nectarchain/user_scripts/vmarandon/FitSPEFFMultiProcess.py @@ -0,0 +1,410 @@ +try: + import sys + import numpy as np + import numpy.ma as ma + from scipy.stats import poisson, norm + from DataUtils import GetLongRunTimeEdges, CountEventTypes, GetTotalEventCounts + from Utils import GetCamera, SignalIntegration, IntegrationMethod, ConvertTitleToName, CustomFormatter, GetDefaultDataPath + from ctapipe.containers import EventType + from FileHandler import DataReader, GetNectarCamEvents + from Stats import CameraSampleStats, CameraStats + from tqdm import tqdm + from FitUtils import JPT2FitFunction + from iminuit import Minuit + import time + import copy + from matplotlib import pyplot as plt + from ctapipe.visualization import CameraDisplay + #from multiprocessing import Pool + import pickle + #from multiprocessing.dummy import Pool as ThreadPool + from multiprocessing import Pool + import argparse + + from IPython import embed +except ImportError as e: + print(e) + raise SystemExit + + +# %matplotlib tk +#ma.array(t0max[0],mask=css.get_lowcount_mask()[0,:,0]) +class FitResult(): + def __init__(self,res,params): + self.minuit_result = res + self.fit_parameters = params + +def split_array(a, n): + k, m = divmod(len(a), n) + return (a[i*k+min(i, m):(i+1)*k+min(i+1, m)] for i in range(n)) + +def save_data(datas, fileName): + """Save info in a file using pickle""" + with open(fileName,'wb') as file: + pickle.dump(datas,file) + +def FitPart(pixels,charges,sigpeds,peds): + results = list() + for p,c,s,ped in tqdm(zip(pixels,charges,sigpeds,peds)): + res = FitPixel(c,s,ped) + if res is not None: + results.append( (p,) + res ) + return results + + +def FitPixel(hicharges,sigped=None,ped=None): + start = time.time() + pix_datas = hicharges + bins = np.linspace( np.min(pix_datas)-0.5,np.max(pix_datas)+0.5,int(np.max(pix_datas)-np.min(pix_datas)+2)) + vals,x_edges = np.histogram(pix_datas,bins=bins) + + if len(bins)<=100: + #print(f"Pixel: {pix} --> Too few bins... bad pixel ? --> SKIP !") + return #continue + + jpt2 = JPT2FitFunction(x_edges,vals) + amplitude = float(np.sum(vals)) + if ped is None: + hiPed = x_edges[ np.argmax(vals) ] + 0.5 + else: + hiPed = ped + + if sigped is None: + sigmaHiPed = 14. + else: + sigmaHiPed = sigped + + ns = 1. + meanillu = 1. + gain = 65. + sigmaGain = 0.45*gain + start_parameters = [amplitude,hiPed,sigmaHiPed,ns,meanillu,gain,sigmaGain] + + m = Minuit( jpt2.Minus2LogLikelihood, start_parameters, name=("Norm","Ped","SigmaPed","Ns","Illu","Gain","SigmaGain") ) + + m.errors["Norm"] = 0.1*amplitude + m.errors["Ped"] = 1.* sigmaHiPed + m.errors["SigmaPed"] = 0.2*sigmaHiPed + m.errors["Ns"] = 0.1 + m.errors["Illu"] = 0.3 + m.errors["Gain"] = 0.3*gain + m.errors["SigmaGain"] = 0.3*sigmaGain + + m.limits['Norm'] = (0.,None) + m.limits['Illu'] = (0.,8.) + m.limits['Gain'] = (10.,200.) + m.limits['SigmaGain'] = (10.,None) + m.limits['Ns'] = (0.,None) + m.limits['SigmaPed'] = (0.,None) + + try: + min_result = m.migrad() + fit_params = m.params + + end = time.time() + fitTime = end-start + #pixelFitTime[pix] = end-start + #pixelFitResult[pix] = FitResult( copy.deepcopy(min_result), copy.deepcopy(fit_params) ) + return FitResult( copy.deepcopy(min_result), copy.deepcopy(fit_params) ), fitTime + + except ValueError as err: + #print(f"Pixel: {pix} Error: {err}") + pass + + +def DoSPEFFFit(arglist): + + p = argparse.ArgumentParser(description='Perform SPE Fit for FF data', + epilog='examples:\n' + '\t python %(prog)s --run 3750 \n', + formatter_class=CustomFormatter) + + p.add_argument("--run", dest='run', type=int, help="Run number") + p.add_argument("--data-path",dest='dataPath',type=str,default=GetDefaultDataPath(),help="Path to the rawdata directory. The program will recursively search in all directory for matching rawdata") + p.add_argument("--njobs", dest='njobs', type=int, default=1, help="Number of CPU to use") + + args = p.parse_args(arglist) + + if args.run is None: + p.print_help() + return -1 + + run = args.run + path = args.dataPath + + # path = '/Users/vm273425/Programs/NectarCAM/data/' + + data = DataReader(run,path) + ok = data.Connect("trigger","r0","mon") + if not ok: + data = GetNectarCamEvents(run,path,applycalib=False) + + nEvents = GetTotalEventCounts(run,path) + + css = CameraSampleStats() + + + counter = 0 + nUsedEntries = 10000000 + for evt in tqdm(data,total=min(nEvents,nUsedEntries)): + css.add( evt.r0.tel[0].waveform, validmask = ~evt.mon.tel[0].pixel_status.hardware_failing_pixels) + counter += 1 + if counter>=nUsedEntries: + break + + mean_waveform = css.mean + + fig1 = plt.figure() + pixs = [777,747,320,380,123,1727,427,74] + for p in pixs: + plt.plot( mean_waveform[0,p,:], label=f'{p}') + plt.grid() + plt.legend() + + + + fig2 = plt.figure() + #pixs = [777,747,320,380,123,1727,427,74] + pixs = [923,483,1573,1751,1491,482,720] + for p in pixs: + plt.plot( mean_waveform[0,p,:], label=f'{p}') + plt.grid() + plt.legend() + + + t0max = np.argmax(mean_waveform,axis=2) + + fig3 = plt.figure() + cam_tom_hg = CameraDisplay(geometry=GetCamera(),image=ma.array(t0max[0],mask=css.get_lowcount_mask()[0,:,0]),cmap='coolwarm',title='High-Gain Mean ToM',allow_pick=True) + cam_tom_hg.highlight_pixels(range(1855),color='grey') + cam_tom_hg.add_colorbar() + #cam_tom_hg.show() + + + fig4 = plt.figure() + histdata = t0max[0] + bins = np.linspace(-0.05,50.05,502) + good_histdata = histdata[ histdata>0. ] + _ = plt.hist(good_histdata,bins=bins) + plt.xlim(np.min(good_histdata)-1,np.max(good_histdata)+1) + plt.title("ToM Average High Gain Waveform") + plt.grid() + + + data = DataReader(run,path) + ok = data.Connect("trigger","r0","mon") + if not ok: + data = GetNectarCamEvents(run,path,applycalib=False) + + nEvents = GetTotalEventCounts(run,path) + + camera = GetCamera() + + + ped_stats = CameraStats() + hicharges = list() + counter = 0 + + #data.rewind() + + counter = 0 + for evt in tqdm(data,total=nEvents): + charges, times = SignalIntegration( evt.r0.tel[0].waveform, exclusion_mask=evt.mon.tel[0].pixel_status.hardware_failing_pixels, peakpositions=t0max,method=IntegrationMethod.USERPEAK,left_bound=5,right_bound=11,camera=camera) + + ped = np.sum( evt.r0.tel[0].waveform[:,:,:16], axis=2 ) + ped_stats.add( ped, validmask=~evt.mon.tel[0].pixel_status.hardware_failing_pixels ) + + hicharges.append( charges[0] ) + counter += 1 +# if counter >=1000: +# break + + hicharges = np.array(hicharges) + + + + fig5 = plt.figure() + cam_pedw_hg = CameraDisplay(geometry=GetCamera(),image=ma.array(ped_stats.stddev[0],mask=ped_stats.get_lowcount_mask()[0]),cmap='coolwarm',title='High-Gain Mean Ped Width',allow_pick=True) + cam_pedw_hg.highlight_pixels(range(1855),color='grey') + cam_pedw_hg.add_colorbar() + #cam_pedw_hg.show() + + + + epedwidth = ped_stats.stddev[0] + eped = ped_stats.mean[0] + badpedwidth = ped_stats.get_lowcount_mask()[0] + + pixelFitResult = dict() + pixelFitTime = dict() + + use_multiprocess = args.njobs>1 + simple_multiprocess = False + # HiCharges = hicharges + # BadPedWidth = badpedwidth + # EPedWidth = epedwidth + + start_fit = time.time() + print(f"Starting the fit at : {start_fit}") + if use_multiprocess and simple_multiprocess: + print("Using the fit by pixel method") + pixlist = [pix for pix in range(1855)] + + with Pool(args.njobs) as p: + datas = list() + sigpeds = list() + peds = list() + for pix in pixlist: + datas.append( hicharges[:,pix] ) + sigpeds.append( 14. if badpedwidth[pix] else epedwidth[pix] ) + peds.append( None if badpedwidth[pix] else eped[pix] ) + + print(len(datas)) + results = p.starmap(FitPixel,zip(datas,sigpeds,peds)) + + p.close() + p.join() + print(len(results)) + + for pix in pixlist: + res = results[pix] + if res is not None: + pixelFitResult[pix] = res[0] + pixelFitTime[pix] = res[1] + elif use_multiprocess and not simple_multiprocess: + print("Using the fit by part method") + nPart = args.njobs + + npix2use = 1855 + part_pixels = list( split_array( [pix for pix in range(npix2use)], nPart ) ) + + part_charges = list( split_array( [hicharges[:,pix] for pix in range(npix2use) ], nPart ) ) + part_sigpeds = list( split_array( [ 14. if badpedwidth[pix] else epedwidth[pix] for pix in range(npix2use)], nPart ) ) + part_peds = list( split_array( [ None if badpedwidth[pix] else eped[pix] for pix in range(npix2use)], nPart ) ) + + #embed() + #for part in part_pixels: + # pcharges + # for pix in part: + with Pool(nPart) as p: + results = p.starmap(FitPart,zip(part_pixels,part_charges,part_sigpeds,part_peds)) + for rs in results: + for res in rs: + if res is not None: + cur_pixel = res[0] + pixelFitResult[cur_pixel] = res[1] + pixelFitTime[cur_pixel] = res[2] + + + + + else: + print("Using the standard approach") + for pix in tqdm(range(1855)): + sigped = 14. if badpedwidth[pix] else epedwidth[pix] + ped = None if badpedwidth[pix] else eped[pix] + res = FitPixel(hicharges[:,pix],sigped,ped ) + if res is not None: + pixelFitResult[pix] = res[0] + pixelFitTime[pix] = res[1] + #if pix>=100 : + # break + + end_fit = time.time() + + print(f"Finishing the fit at : {end_fit}") + print(f"Time to do it : {end_fit-start_fit} s") + + + print(len(pixelFitResult)) + + gains = list() + for p, f in pixelFitResult.items(): + fit_gain = f.fit_parameters['Gain'].value + gains.append(fit_gain) + + + cam_gain_values = np.zeros(1855) + for p, f in pixelFitResult.items(): + fit_gain = f.fit_parameters['Gain'].value + cam_gain_values[p] = fit_gain + + + fig6 = plt.figure(figsize=(8,8)) + cam_gain = CameraDisplay(geometry=GetCamera(),image=ma.array(cam_gain_values,mask=cam_gain_values<40.),cmap='turbo',title=f'High-Gain Gain (VIM Edition)\nRun: {run}',show_frame=False) + cam_gain.highlight_pixels(range(1855),color='grey') + cam_gain.add_colorbar() + cam_gain.show() + #fig.savefig(f'run_{run}_gain_camera.png') + fig6.savefig(f'run_{run}_FittedGain_VIMEdition.png') + + fig7 = plt.figure(figsize=(8,8)) + mean_gain = np.mean(gains) + std_gain = np.std(gains) + _ = plt.hist(gains,bins=60,label=f'$\mu:$ {mean_gain:.2f}\n$\sigma:$ {std_gain:.2f}') + plt.title(f"Fitted Gain Distribution (Run: {run})") + plt.xlabel('Gain (ADC per pe)') + plt.grid() + plt.legend() + #fig.savefig(f'run_{run}_gain_distribution.png') + fig7.savefig(f'run_{run}_FittedGainDistribution_VIMEdition.png') + + cam_illu_values = np.zeros(1855) + for p, f in pixelFitResult.items(): + fit_illu = f.fit_parameters['Illu'].value + cam_illu_values[p] = fit_illu + + fig8 = plt.figure(figsize=(8,8)) + cam_illu = CameraDisplay(geometry=GetCamera(),image=ma.array(cam_illu_values,mask=cam_gain_values<40.),cmap='turbo',title=f'Mean Illumination (VIM Edition)\nRun: {run}',show_frame=False) + cam_illu.highlight_pixels(range(1855),color='grey') + cam_illu.add_colorbar() + cam_illu.show() + fig8.savefig(f'run_{run}_FittedIllumination_VIMEdition.png') + + cam_ns_values = np.zeros(1855) + for p, f in pixelFitResult.items(): + fit_ns = f.fit_parameters['Ns'].value + cam_ns_values[p] = fit_ns + + fig9 = plt.figure(figsize=(8,8)) + cam_ns = CameraDisplay(geometry=GetCamera(),image=ma.array(cam_ns_values,mask=cam_gain_values<40.),cmap='turbo',title=f'Ns Parameter (VIM Edition)\nRun: {run}',show_frame=False) + cam_ns.highlight_pixels(range(1855),color='grey') + cam_ns.add_colorbar() + cam_ns.show() + + cam_pedw_values = np.zeros(1855) + for p, f in pixelFitResult.items(): + fit_pedw = f.fit_parameters['SigmaPed'].value + cam_pedw_values[p] = fit_pedw + + fig9 = plt.figure(figsize=(8,8)) + cam_pedw = CameraDisplay(geometry=GetCamera(),image=ma.array(cam_pedw_values,mask=cam_gain_values<40.),cmap='turbo',title=f'Pedestal Width (VIM Edition)\nRun: {run}',show_frame=False) + cam_pedw.highlight_pixels(range(1855),color='grey') + cam_pedw.add_colorbar() + cam_pedw.show() + fig9.savefig(f'run_{run}_FittedPedestalWidth_VIMEdition.png') + + + ## Save results : + outdataname = f'run_{run}_FittedSPEResult.dat' + + with open(outdataname,'w') as outfile: + outfile.write(f"pix Ped SigmaPed Ns Illu Gain SigmaGain\n") + for p, f in pixelFitResult.items(): + if f.minuit_result.valid: + outfile.write(f"{p} {f.fit_parameters['Ped'].value} {f.fit_parameters['SigmaPed'].value} {f.fit_parameters['Ns'].value} {f.fit_parameters['Illu'].value} {f.fit_parameters['Gain'].value} {f.fit_parameters['SigmaGain'].value}\n") + + # m.errors["Norm"] = 0.1*amplitude + # m.errors["Ped"] = 1.* sigmaHiPed + # m.errors["SigmaPed"] = 0.2*sigmaHiPed + # m.errors["Ns"] = 0.1 + # m.errors["Illu"] = 0.3 + # m.errors["Gain"] = 0.3*gain + # m.errors["SigmaGain"] = 0.3*sigmaGain + + plt.show() + + +if __name__ == '__main__': + DoSPEFFFit(sys.argv[1:]) + diff --git a/src/nectarchain/user_scripts/vmarandon/FitUtils.py b/src/nectarchain/user_scripts/vmarandon/FitUtils.py new file mode 100644 index 00000000..199203d6 --- /dev/null +++ b/src/nectarchain/user_scripts/vmarandon/FitUtils.py @@ -0,0 +1,95 @@ +try: + import sys + import numpy as np + from scipy.stats import poisson, norm + + from IPython import embed + +except ImportError as e: + print(e) + raise SystemExit + + +class GaussHistoFitFunction(): + def __init__(self,bin_edges,bin_contents,integrate=True): # Assume that the bin centers are "integers" and that they are spaced at 1.... refinement will be for later + self.xcentre = 0.5 * (bin_edges[1:] + bin_edges[:-1]) + self.xedge = bin_edges + self.y = bin_contents + self.use_integration = integrate + + def ExpectedValues(self,normalisation,mu,sigma): + y = normalisation*np.exp(-0.5*np.power((self.xcentre-mu)/sigma,2.))/(np.sqrt(2.*np.pi)*sigma) + return y + + + def IntegratedExpectedValues(self,normalisation,mu,sigma): + yup = norm.cdf(self.xedge[1:],loc=mu,scale=sigma) + ylow = norm.cdf(self.xedge[:-1],loc=mu,scale=sigma) + return normalisation*(yup-ylow) + + + def Minus2LogLikelihood(self,parameters): + normalisation = parameters[0] + mu = parameters[1] + sigma = parameters[2] + if self.use_integration: + y_fit = self.IntegratedExpectedValues(normalisation,mu,sigma) + else: + y_fit = self.ExpectedValues(normalisation,mu,sigma) + #print(y_fit) + log_likes = poisson.logpmf(self.y,y_fit) + total = np.sum(log_likes) + return -2.*total + + +class JPT2FitFunction(): + def __init__(self,bin_edges,bin_contents): + self.xcentre = 0.5 * (bin_edges[1:] + bin_edges[:-1]) + self.xedge = bin_edges + self.y = bin_contents.astype(float) + + def Prediction(self,parameters): + amplitude = parameters[0] + hiPed = parameters[1] + sigmaHiPed = parameters[2] + ns = parameters[3] + meanillu = parameters[4] + gain = parameters[5] + sigmaGain = parameters[6] + + ## Compute the maximum illumination to consider + nMuUsed = int(5.*np.sqrt(meanillu+1)) + 2 + + ## Precompute some parameters that will be needed for the full iteration + imu = np.arange(nMuUsed) + sig2 = imu * (sigmaGain * sigmaGain) + sigmaHiPed*sigmaHiPed + pos_x = imu * gain + hiPed + poiG = poisson.pmf(imu,meanillu) / np.sqrt(2.*np.pi*sig2) + poiG[1:] *= ns + + sig2x2 = 2.*sig2 + ## Loop on pixels: + + + predicted_entries = np.zeros_like( self.y ) + + #embed() + for i in range(self.xcentre.shape[0]): + predicted_entries[i] = np.sum( poiG * np.exp( -np.power( self.xcentre[i] - pos_x,2.)/sig2x2 ) ) + + predicted_entries *= amplitude + return predicted_entries + + + def Minus2LogLikelihood(self,parameters): + #print(f'Minuit2 Called {parameters}') + y_pred = self.Prediction(parameters) + + LogLikeBins = poisson.logpmf(self.y,y_pred) + LogLike = np.sum( LogLikeBins ) + #if np.isnan(LogLike): + # print("--> Nan") + Minus2LogLike = -2.*LogLike + return Minus2LogLike + + diff --git a/src/nectarchain/user_scripts/vmarandon/ShowHardwareProblem.py b/src/nectarchain/user_scripts/vmarandon/ShowHardwareProblem.py new file mode 100644 index 00000000..a902519c --- /dev/null +++ b/src/nectarchain/user_scripts/vmarandon/ShowHardwareProblem.py @@ -0,0 +1,166 @@ +try: + import sys + + import numpy as np + import numpy.ma as ma + + import argparse + + import matplotlib + from matplotlib import pyplot as plt + + from FileHandler import DataReader, GetNectarCamEvents + from DataUtils import GetTotalEventCounts + from Utils import GetEventTypeFromString, GetCamera, CustomFormatter, GetDefaultDataPath + from FitUtils import GaussHistoFitFunction + from Stats import CameraSampleStats, CameraStats + from CalibrationData import CalibrationCameraDisplay + + from scipy.optimize import curve_fit + from scipy.stats import norm + from iminuit import Minuit + + from ctapipe.visualization import CameraDisplay + + from tqdm import tqdm + + from IPython import embed + +except ImportError as e: + print(e) + raise SystemExit + +class HardwareProblemStatsMaker: + ''' + Class to compute stats on possible hardware problems + ''' + def __init__(self,telid = 0): + self.telid = telid + self.missing = None + +def ShowHardwareProblem(arglist): + + p = argparse.ArgumentParser(description='Show stats on the hardware_failing_pixels flag', + epilog='examples:\n' + '\t python %(prog)s --run 123456 \n', + formatter_class=CustomFormatter) + + p.add_argument("--run", dest='run', type=int, help="Run number to be converted") + p.add_argument("--data-path",dest='dataPath',type=str,default=None,help="Path to the rawdata directory. The program will recursively search in all directory for matching rawdata") + p.add_argument("--telid",dest="telId",type=int,default=0,help="Telescope id for which we show the data distributions") + p.add_argument("--list",dest="listEvent",action='store_true',help="List pixels for which there is hardware problems") + p.add_argument("--savefig",dest="savefig",action='store_true',help='Save figure') + p.add_argument("--debug",dest="debug",action='store_true',help='Debug mode (only 10000 events processed)') + + args = p.parse_args(arglist) + + if args.run is None: + p.print_help() + return -1 + + if args.dataPath is None: + path = GetDefaultDataPath() + else: + path = args.dataPath + + data = DataReader(args.run,path) + ok = data.Connect( "mon","trigger" ) + #evt.mon.tel[self.telid].pixel_status.hardware_failing_pixels + if not ok: + data = GetNectarCamEvents(args.run,path,applycalib=False) + + dataevents = GetTotalEventCounts(args.run,path) + + camstats = CameraStats() + + sum_failing = None + counter = 0 + try: + for evt in tqdm(data,total=dataevents): + hardware_failing = evt.mon.tel[args.telId].pixel_status.hardware_failing_pixels + camstats.add( hardware_failing ) + + if sum_failing is not None: + sum_failing += hardware_failing + else: + sum_failing = hardware_failing.copy().astype(int) + + + if args.debug: + counter += 1 + if counter > 10000: + break + except Exception as e: + print(e) + + lw = 0.2 + + fig_counts, axs_counts = plt.subplots(nrows=1,ncols=2, figsize=(12,6)) + + cam_count_hg = CameraDisplay(geometry=GetCamera(), cmap='turbo', image=ma.array(camstats.mean[0]*camstats.count[0],mask=camstats.mean[0]==0.), title=f'HG Hardware Failing Counts\nrun {args.run}',ax=axs_counts[0],show_frame=False,allow_pick=True) + cam_count_hg.highlight_pixels(range(1855),linewidth=lw,color='grey') + cam_count_hg.add_colorbar() + + cam_count_lg = CameraDisplay(geometry=GetCamera(), cmap='turbo', image=ma.array(camstats.mean[1]*camstats.count[1],mask=camstats.mean[1]==0.), title=f'LG Hardware Failing Counts\nrun {args.run}',ax=axs_counts[1],show_frame=False,allow_pick=True) + cam_count_lg.highlight_pixels(range(1855),linewidth=lw,color='grey') + cam_count_lg.add_colorbar() + + if args.savefig: + figname = f'run_{args.run}_Hardware_Failing_Pixels_Counter.png' + + fig_counts.show() + + + fig_counts2, axs_counts2 = plt.subplots(nrows=1,ncols=2, figsize=(12,6)) + + cam_count_hg2 = CameraDisplay(geometry=GetCamera(), cmap='turbo', image=ma.array(sum_failing[0],mask=sum_failing[0]==0.), title=f'HG Hardware Failing Counts 2\nrun {args.run}',ax=axs_counts2[0],show_frame=False,allow_pick=True) + cam_count_hg2.highlight_pixels(range(1855),linewidth=lw,color='grey') + cam_count_hg2.add_colorbar() + + cam_count_lg2 = CameraDisplay(geometry=GetCamera(), cmap='turbo', image=ma.array(sum_failing[1],mask=sum_failing[1]==0.), title=f'LG Hardware Failing Counts 2\nrun {args.run}',ax=axs_counts2[1],show_frame=False,allow_pick=True) + cam_count_lg2.highlight_pixels(range(1855),linewidth=lw,color='grey') + cam_count_lg2.add_colorbar() + + fig_counts2.show() + + + + fig_frac, axs_frac = plt.subplots(nrows=1,ncols=2, figsize=(12,6)) + + cam_frac_hg = CameraDisplay(geometry=GetCamera(), cmap='turbo', image=ma.array(camstats.mean[0]*100,mask=camstats.mean[0]==0.), title=f'HG Hardware Failing Fraction (%)\nrun {args.run}',ax=axs_frac[0],show_frame=False,allow_pick=True) + cam_frac_hg.add_colorbar() + cam_frac_hg.highlight_pixels(range(1855),linewidth=lw,color='grey') + cam_frac_hg.colorbar.set_label('%') + + cam_frac_lg = CameraDisplay(geometry=GetCamera(), cmap='turbo', image=ma.array(camstats.mean[1]*100,mask=camstats.mean[1]==0.), title=f'LG Hardware Failing Fraction (%)\nrun {args.run}',ax=axs_frac[1],show_frame=False,allow_pick=True) + cam_frac_lg.add_colorbar() + cam_frac_lg.highlight_pixels(range(1855),linewidth=lw,color='grey') + cam_frac_lg.colorbar.set_label('%') + + if args.savefig: + figname = f'run_{args.run}_Hardware_Failing_Pixels_Fraction.png' + + fig_frac.show() + + + #embed() + + ## To go further one could use this information : + #module_status[ evt.nectarcam.tel[0].svc.module_ids ] = evt.nectarcam.tel[0].evt.module_status + # to know the module that was expected to be there during acquisition in order to avoid false positive. + # There is a similar info for pixels + #pixel_status = np.zeros(N_PIXELS) + #pixel_status[self.camera_config.expected_pixels_id] = event.pixel_status + #status_container.hardware_failing_pixels[:] = pixel_status == 0 + + + + #plt.show() + input("Press Enter to quit\n") + + embed() + + + +if __name__ == "__main__": + ShowHardwareProblem(sys.argv[1:]) diff --git a/src/nectarchain/user_scripts/vmarandon/Utils.py b/src/nectarchain/user_scripts/vmarandon/Utils.py index a433ddc6..fe4af69c 100644 --- a/src/nectarchain/user_scripts/vmarandon/Utils.py +++ b/src/nectarchain/user_scripts/vmarandon/Utils.py @@ -6,6 +6,9 @@ import datetime import re import enum + import pickle + import lz4.frame + from ctapipe.io import EventSource from ctapipe.instrument import CameraGeometry @@ -37,6 +40,26 @@ class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawDescriptionHelpFormatter): pass +def ConvertTitleToName(title: str): + title = title.replace(' ','_') + title = title.replace("\n","_") + title = title.replace("(","_") + title = title.replace(")","_") + title = title.replace(":","_") + title = title.replace(",","_") + title = title.replace("=","_") + title = title.replace("_-_","_") + title = title.replace(">","Above") + title = title.replace("<","Below") + #title = title.replace("0.","0d") + title = title.replace(".","d") + + + ## do a bit of cleanup of the _ + while title.find('__') != -1: + title = title.replace("__","_") + return title + def GetDefaultDataPath(): return os.environ.get( 'NECTARCAMDATA' , '/Users/vm273425/Programs/NectarCAM/data') @@ -145,11 +168,12 @@ def clean_waveform(wvf, use_nan = True): class IntegrationMethod(enum.Enum): PEAKSEARCH = enum.auto() NNSEARCH = enum.auto() + USERPEAK = enum.auto() -def SignalIntegration(waveform,exclusion_mask=None, method = IntegrationMethod.PEAKSEARCH,left_bound=5,right_bound=7,camera=None): +def SignalIntegration(waveform,exclusion_mask=None, method = IntegrationMethod.PEAKSEARCH,left_bound=5,right_bound=7,camera=None,peakpositions=None): wvf = waveform.astype(float) if exclusion_mask is not None: wvf[ exclusion_mask ] = 0. @@ -160,6 +184,9 @@ def SignalIntegration(waveform,exclusion_mask=None, method = IntegrationMethod.P if camera is None: camera = GetCamera() charge, time = NeighborPeakIntegration(waveform, camera.neighbor_matrix, left_bound=5, right_bound=7) + elif method == IntegrationMethod.USERPEAK: + #print("HERE") + charge, time = UserPeakIntegration(waveform,peakpositions=peakpositions,left_bound=left_bound,right_bound=right_bound) else: print("I Don't know about this method !") @@ -229,6 +256,37 @@ def NeighborPeakIntegration(waveform, neighbor_matrix, left_bound=5, right_bound return integrated_signal, signal_timeslice +@njit(parallel=True) +def UserPeakIntegration(waveform, peakpositions, left_bound=5, right_bound=7): + + wvf_shape = waveform.shape + n_channel = wvf_shape[0] + n_pixels = wvf_shape[1] + n_samples = wvf_shape[2] + + integrated_signal = np.zeros( (n_channel,n_pixels) ) + signal_timeslice = np.zeros( (n_channel,n_pixels) ) + + integ_window = left_bound+right_bound + + for chan in prange(2): + chan_wvf = waveform[chan,:,:] + for pix in range(n_pixels): + trace = chan_wvf[pix,:] + peak_pos = peakpositions[chan,pix] + if (peak_pos-left_bound) < 0: + lo_bound = 0 + up_bound = integ_window + elif (peak_pos + right_bound) > n_samples: + lo_bound = n_samples-integ_window + up_bound = n_samples + else: + lo_bound = peak_pos - left_bound + up_bound = peak_pos + right_bound + integrated_signal[chan,pix] = np.sum( trace[lo_bound:up_bound] ) + signal_timeslice[chan,pix] = peak_pos + + return integrated_signal, signal_timeslice def getPixelT0Spline(waveform): @@ -368,3 +426,17 @@ def GetEventTypeFromString(event_str): print(f"WARNING> Don't know about the event type [{event_str}]") evt_type = EventType.UNKNOWN return evt_type + + +def save_simple_data(datas, fileName): + """Save info in a file using pickle""" + with lz4.frame.open( fileName, 'wb' ) as f : + pickle.dump(datas,f) + + #with open(fileName,'wb') as file: + # pickle.dump(datas,file) + +def read_simple_data(fileName): + """Read info from a file using pickle""" + with lz4.frame.open( fileName, 'rb' ) as f : + return pickle.load(f)