Coverage for lst_auto_rta/Plot_event_DL1.py: 0%

1#!/usr/bin/env python3 

2import argparse 

3import os 

4from os import listdir 

5from os.path import isfile, join 

6 

7import matplotlib.pyplot as plt 

8import numpy as np 

9import tables 

10from matplotlib.backends.backend_pdf import PdfPages 

11 

12parser = argparse.ArgumentParser( 

13 description="Automatic Script for the DL1 check", formatter_class=argparse.ArgumentDefaultsHelpFormatter 

14) 

15parser.add_argument("-d", "--directory", default="/fefs/onsite/pipeline/rta/data/", help="Directory for data") 

16parser.add_argument("-da", "--date", default="20230705", help="Date of the run to check") 

17parser.add_argument("-r", "--run-id", default="13600", help="run id to check") 

18parser.add_argument("-add", "--add-string", default="reco/", help="add a string to the path") 

19parser.add_argument("-o", "--output", default="", help="output path") 

20 

21args = parser.parse_args() 

22config = vars(args) 

23 

24font = {"family": "normal", "weight": "bold", "size": 22} 

25 

26plt.rc("font", **font) 

27 

28 

29def plotImage(outputPlot, tabX, tabY, signal, imageIndex, eventId, eventType, tabFixPixelOrder=None, id_=0): 

30 strInfo = f"image {imageIndex}, id = {eventId}, type = {eventType},\n Signal mean = {signal.mean()} pe, min = {signal.min()}, max = {signal.max()}, std = {signal.std()}" 

31 print(strInfo) 

32 fig = plt.figure(figsize=(16, 10)) 

33 fig.patch.set_alpha(1.0) 

34 if tabFixPixelOrder is not None: 

35 plt.scatter(tabX[tabFixPixelOrder], tabY[tabFixPixelOrder], c=signal, s=120) 

36 else: 

37 plt.scatter(tabX, tabY, c=signal, s=120) 

38 

39 plt.axis("equal") 

40 plt.xlabel("x") 

41 plt.ylabel("y") 

42 plt.colorbar() 

43 plt.text(0.05, 0.95, strInfo, transform=fig.transFigure, size=14) 

44 plt.show() 

45 plt.savefig(outputPlot + "/plotImage" + str(id_) + ".png") 

46 # pdf.savefig() # saves the current figure into a pdf page 

47 plt.close() 

48 

49 

50def processAllFile(outputPlot, inputFile, firstImageIndex=0, lastImageIndex=10000, indexStep=1, id_=0): 

51 hfile = tables.open_file(inputFile, "r") 

52 

53 tabX = hfile.root.configuration.instrument.telescope.camera.geometry_LSTCam.col("pix_x") 

54 tabY = hfile.root.configuration.instrument.telescope.camera.geometry_LSTCam.col("pix_y") 

55 

56 try: 

57 tabPixelOrder = hfile.root.configuration.instrument.telescope.camera.pixel_order.tel_001.read() 

58 except: 

59 tabPixelOrder = None 

60 

61 # tabImages = hfile.root.dl1.event.telescope.images.tel_001.col("image") 

62 # tabEventId = hfile.root.dl1.event.subarray.trigger.col("event_id") 

63 # tabEventType = hfile.root.dl1.event.subarray.trigger.col("event_type") 

64 

65 tabImages = hfile.root.dl1.event.telescope.images.tel_001.col("image") 

66 tabTime = hfile.root.dl1.event.telescope.images.tel_001.col("peak_time") 

67 tabEventId = hfile.root.dl1.event.telescope.parameters.tel_001.col("event_id") 

68 tabEventType = hfile.root.dl1.event.telescope.parameters.tel_001.col("event_type") 

69 tabIntensity = hfile.root.dl1.event.telescope.parameters.tel_001.col("hillas_intensity") 

70 # id_ = 0 

71 # for i in range(firstImageIndex, len(tabImages), indexStep): 

72 for i in range(firstImageIndex, lastImageIndex, indexStep): 

73 if tabEventType[i] == 1 and tabIntensity[i] > 11000 and tabIntensity[i] < 12000: 

74 plotImage(outputPlot, tabX, tabY, tabImages[i], i, tabEventId[i], tabEventType[i], tabPixelOrder, id_) 

75 id_ = id_ + 1 

76 return id_ 

77 

78 

79mypath = config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] + "/dl1/" 

80mypath_rundir = config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] 

81print("Path is here ", mypath) 

82if not os.path.exists(mypath_rundir + "plots"): 

83 os.mkdir(mypath_rundir + "plots") 

84filename = [f for f in listdir(mypath) if isfile(join(mypath, f)) if "dl1" in f if not "v06" in f] 

85print("Plotting image from these files", filename) 

86index = 0 

87for i in filename: 

88 print(mypath + i) 

89 print("output in ", mypath_rundir + "plots") 

90 index = processAllFile(mypath_rundir + "plots", mypath + i, id_=index)