From 426b7cb01057cef1f216b1b63eb203b5df007f71 Mon Sep 17 00:00:00 2001 From: Karagiannis Chrysovalantis <102170668+ckaragian@users.noreply.github.com> Date: Fri, 16 Feb 2024 19:25:29 +0200 Subject: [PATCH] Update to Track Length in the water & Energy Reconstruction (#248) * Adding tools for the track length in the water tank and the muon energy reconstruction * Adding the config files for the EnergyReco toolchains * Correcting some mistakes in README.md files * Fixed units of MC track length in mrd in FindTrackLengthInWater tool and converted reconstructed track length in the water and energy in [m] and [GeV] before passing them to the ANNIEEvent store. * Update PlotsTrackLengthAndEnergy.cpp Fix memory leaks * Fixing minor bug in PlotsTrackLengthAndEnergy.cpp script --------- Co-authored-by: marc1uk --- .../BDTMuonEnergyPredict.py | 177 ++++++ UserTools/BDTMuonEnergyPredict/README.md | 40 ++ .../BDTMuonEnergyTrain_Test.py | 280 +++++++++ UserTools/BDTMuonEnergyTrain_Test/README.md | 44 ++ .../DNNTrackLengthPredict.py | 201 +++++++ UserTools/DNNTrackLengthPredict/README.md | 45 ++ .../DNNTrackLengthTrain_Test.py | 335 +++++++++++ UserTools/DNNTrackLengthTrain_Test/README.md | 49 ++ UserTools/Factory/Factory.cpp | 1 + .../FindTrackLengthInWater.cpp | 557 +++++++++++------- .../FindTrackLengthInWater.h | 32 +- UserTools/FindTrackLengthInWater/README.md | 71 ++- .../PlotsTrackLengthAndEnergy.cpp | 121 ++++ .../PlotsTrackLengthAndEnergy.h | 45 ++ UserTools/PlotsTrackLengthAndEnergy/README.md | 23 + UserTools/Unity.h | 1 + configfiles/EnergyReco/DummyToolConfig | 3 - configfiles/EnergyReco/Erecoconfig | 6 - .../Predict/BDTMuonEnergyPredictConfig | 9 + .../EnergyReco/Predict/ClusterFinderConfig | 12 + .../Predict/DNNTrackLengthPredictConfig | 13 + .../EnergyReco/Predict/DigitBuilderConfig | 11 + .../EnergyReco/Predict/EventSelectorConfig | 15 + .../EnergyReco/Predict/FindMrdTracksConfig | 12 + .../Predict/FindTrackLengthInWaterConfig | 5 + .../EnergyReco/Predict/HitCleanerConfig | 20 + configfiles/EnergyReco/Predict/LAPPDIDs.txt | 5 + .../EnergyReco/Predict/LoadWCSimConfig | 15 + .../EnergyReco/Predict/LoadWCSimLAPPDConfig | 10 + .../Predict/MCParticlePropertiesConfig | 1 + .../Predict/MCRecoEventLoaderConfig | 12 + .../Predict/PlotsTrackLengthAndEnergyConfig | 1 + .../EnergyReco/Predict/TimeClusteringConfig | 12 + .../EnergyReco/{ => Predict}/ToolChainConfig | 6 +- configfiles/EnergyReco/Predict/ToolsConfig | 19 + .../Predict/VtxExtendedVertexFinderConfig | 9 + .../EnergyReco/Predict/VtxSeedFineGridConfig | 7 + .../EnergyReco/Predict/VtxSeedGeneratorConfig | 7 + configfiles/EnergyReco/README.md | 57 +- configfiles/EnergyReco/ToolsConfig | 1 - .../Train_Test/BDTMuonEnergyTrainConfig | 15 + .../EnergyReco/Train_Test/ClusterFinderConfig | 12 + .../Train_Test/DNNTrackLengthTrainConfig | 18 + .../EnergyReco/Train_Test/DigitBuilderConfig | 11 + .../EnergyReco/Train_Test/EventSelectorConfig | 15 + .../EnergyReco/Train_Test/FindMrdTracksConfig | 12 + .../Train_Test/FindTrackLengthInWaterConfig | 6 + .../EnergyReco/Train_Test/HitCleanerConfig | 20 + .../EnergyReco/Train_Test/LAPPDIDs.txt | 5 + .../EnergyReco/Train_Test/LoadWCSimConfig | 15 + .../Train_Test/LoadWCSimLAPPDConfig | 10 + .../Train_Test/MCParticlePropertiesConfig | 1 + .../Train_Test/MCRecoEventLoaderConfig | 12 + .../Train_Test/TimeClusteringConfig | 12 + .../EnergyReco/Train_Test/ToolChainConfig | 23 + configfiles/EnergyReco/Train_Test/ToolsConfig | 18 + .../Train_Test/VtxExtendedVertexFinderConfig | 9 + .../Train_Test/VtxSeedFineGridConfig | 7 + .../Train_Test/VtxSeedGeneratorConfig | 6 + 59 files changed, 2264 insertions(+), 253 deletions(-) create mode 100644 UserTools/BDTMuonEnergyPredict/BDTMuonEnergyPredict.py create mode 100644 UserTools/BDTMuonEnergyPredict/README.md create mode 100644 UserTools/BDTMuonEnergyTrain_Test/BDTMuonEnergyTrain_Test.py create mode 100644 UserTools/BDTMuonEnergyTrain_Test/README.md create mode 100644 UserTools/DNNTrackLengthPredict/DNNTrackLengthPredict.py create mode 100644 UserTools/DNNTrackLengthPredict/README.md create mode 100644 UserTools/DNNTrackLengthTrain_Test/DNNTrackLengthTrain_Test.py create mode 100644 UserTools/DNNTrackLengthTrain_Test/README.md create mode 100644 UserTools/PlotsTrackLengthAndEnergy/PlotsTrackLengthAndEnergy.cpp create mode 100644 UserTools/PlotsTrackLengthAndEnergy/PlotsTrackLengthAndEnergy.h create mode 100644 UserTools/PlotsTrackLengthAndEnergy/README.md delete mode 100644 configfiles/EnergyReco/DummyToolConfig delete mode 100644 configfiles/EnergyReco/Erecoconfig create mode 100644 configfiles/EnergyReco/Predict/BDTMuonEnergyPredictConfig create mode 100644 configfiles/EnergyReco/Predict/ClusterFinderConfig create mode 100644 configfiles/EnergyReco/Predict/DNNTrackLengthPredictConfig create mode 100644 configfiles/EnergyReco/Predict/DigitBuilderConfig create mode 100644 configfiles/EnergyReco/Predict/EventSelectorConfig create mode 100644 configfiles/EnergyReco/Predict/FindMrdTracksConfig create mode 100644 configfiles/EnergyReco/Predict/FindTrackLengthInWaterConfig create mode 100644 configfiles/EnergyReco/Predict/HitCleanerConfig create mode 100644 configfiles/EnergyReco/Predict/LAPPDIDs.txt create mode 100644 configfiles/EnergyReco/Predict/LoadWCSimConfig create mode 100644 configfiles/EnergyReco/Predict/LoadWCSimLAPPDConfig create mode 100644 configfiles/EnergyReco/Predict/MCParticlePropertiesConfig create mode 100644 configfiles/EnergyReco/Predict/MCRecoEventLoaderConfig create mode 100644 configfiles/EnergyReco/Predict/PlotsTrackLengthAndEnergyConfig create mode 100644 configfiles/EnergyReco/Predict/TimeClusteringConfig rename configfiles/EnergyReco/{ => Predict}/ToolChainConfig (71%) create mode 100644 configfiles/EnergyReco/Predict/ToolsConfig create mode 100644 configfiles/EnergyReco/Predict/VtxExtendedVertexFinderConfig create mode 100644 configfiles/EnergyReco/Predict/VtxSeedFineGridConfig create mode 100644 configfiles/EnergyReco/Predict/VtxSeedGeneratorConfig delete mode 100644 configfiles/EnergyReco/ToolsConfig create mode 100644 configfiles/EnergyReco/Train_Test/BDTMuonEnergyTrainConfig create mode 100644 configfiles/EnergyReco/Train_Test/ClusterFinderConfig create mode 100644 configfiles/EnergyReco/Train_Test/DNNTrackLengthTrainConfig create mode 100644 configfiles/EnergyReco/Train_Test/DigitBuilderConfig create mode 100644 configfiles/EnergyReco/Train_Test/EventSelectorConfig create mode 100644 configfiles/EnergyReco/Train_Test/FindMrdTracksConfig create mode 100644 configfiles/EnergyReco/Train_Test/FindTrackLengthInWaterConfig create mode 100644 configfiles/EnergyReco/Train_Test/HitCleanerConfig create mode 100644 configfiles/EnergyReco/Train_Test/LAPPDIDs.txt create mode 100644 configfiles/EnergyReco/Train_Test/LoadWCSimConfig create mode 100644 configfiles/EnergyReco/Train_Test/LoadWCSimLAPPDConfig create mode 100644 configfiles/EnergyReco/Train_Test/MCParticlePropertiesConfig create mode 100644 configfiles/EnergyReco/Train_Test/MCRecoEventLoaderConfig create mode 100644 configfiles/EnergyReco/Train_Test/TimeClusteringConfig create mode 100644 configfiles/EnergyReco/Train_Test/ToolChainConfig create mode 100644 configfiles/EnergyReco/Train_Test/ToolsConfig create mode 100644 configfiles/EnergyReco/Train_Test/VtxExtendedVertexFinderConfig create mode 100644 configfiles/EnergyReco/Train_Test/VtxSeedFineGridConfig create mode 100644 configfiles/EnergyReco/Train_Test/VtxSeedGeneratorConfig diff --git a/UserTools/BDTMuonEnergyPredict/BDTMuonEnergyPredict.py b/UserTools/BDTMuonEnergyPredict/BDTMuonEnergyPredict.py new file mode 100644 index 000000000..74f81a19f --- /dev/null +++ b/UserTools/BDTMuonEnergyPredict/BDTMuonEnergyPredict.py @@ -0,0 +1,177 @@ +# BDTMuonEnergy Tool script +# ------------------ +from Tool import * +import sys +import numpy as np +import pandas as pd +import tensorflow as tf +import tempfile +import random +import csv +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from array import array +from sklearn import datasets +from sklearn import metrics +from sklearn import model_selection +from sklearn import preprocessing +import sklearn +from sklearn.utils import shuffle +from sklearn import linear_model, ensemble +from sklearn.metrics import mean_squared_error +import pickle +import ROOT + +class BDTMuonEnergyPredict(Tool): + + # declare member variables here + weightsfilename = std.string("") + + def Initialise(self): + self.m_log.Log(__file__+" Initialising", self.v_debug, self.m_verbosity) + self.m_variables.Print() + self.m_variables.Get("BDTMuonModelFile", self.weightsfilename) + return 1 + + def Execute(self): + self.m_log.Log(__file__+" Executing", self.v_debug, self.m_verbosity) + #Load Data + EnergyRecoBoostStore=cppyy.gbl.BoostStore(True, 2)#define the energy boost store class object to load the variables for the BDT training + EnergyRecoBoostStore = self.m_data.Stores.at("EnergyReco") + #Retrieve the required variables from the Store + EnergyRecoBoostStore.Print(False) + get_ok = EnergyRecoBoostStore.Has("num_pmt_hits") + num_pmt_hits=ctypes.c_int(0) + if not get_ok: + print("BDTMuonEnergyPredict Tool: There is no entry in Energy Reco boost store.") + return 1 + if get_ok: + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry num_pmt_hits: ", get_ok) + print("BDTMuonEnergyPredict Tool: type of num_pmt_hits entry is :",EnergyRecoBoostStore.Type("num_pmt_hits")) + print("BDTMuonEnergyPredict Tool: Getting num_pmt_hits from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("num_pmt_hits",num_pmt_hits) + print("BDTMuonEnergyPredict Tool: Number of pmt hits is: ", num_pmt_hits.value) + ok = EnergyRecoBoostStore.Has("num_lappd_hits") + num_lappd_hits=ctypes.c_int(0) + if ok: + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry num_lappd_hits: ",ok) + print("BDTMuonEnergyPredict Tool: type of num_lappd_hits entry is :",EnergyRecoBoostStore.Type("num_lappd_hits")) + print("BDTMuonEnergyPredict Tool: Getting num_lappd_hits from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("num_lappd_hits",num_lappd_hits) + print("BDTMuonEnergyPredict Tool: Number of lappd hits is: ", num_lappd_hits.value) + ok = EnergyRecoBoostStore.Has("DNNRecoLength") + DNNRecoLength=ctypes.c_double(0) + if ok: + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry DNNRecoLength: ",ok) + print("BDTMuonEnergyPredict Tool: type of DNNRecoLength entry is :",EnergyRecoBoostStore.Type("DNNRecoLength")) + print("BDTMuonEnergyPredict Tool: Getting DNNRecoLength from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("DNNRecoLength",DNNRecoLength) + print("BDTMuonEnergyPredict Tool: The reconstructed track length in the water by the DNN is: ", DNNRecoLength.value) + ok = EnergyRecoBoostStore.Has("recoTrackLengthInMrd") + recoTrackLengthInMrd=ctypes.c_double(0) + if ok: + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry recoTrackLengthInMrd: ",ok) + print("BDTMuonEnergyPredict Tool: type of recoTrackLengthInMrd entry is :",EnergyRecoBoostStore.Type("recoTrackLengthInMrd")) + print("BDTMuonEnergyPredict Tool: Getting recoTrackLengthInMrd from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("recoTrackLengthInMrd",recoTrackLengthInMrd) + print("BDTMuonEnergyPredict Tool: The reconstructed track length in the MRD is: ", recoTrackLengthInMrd.value) + ok = EnergyRecoBoostStore.Has("diffDirAbs") + diffDirAbs=ctypes.c_float(0) + if ok: + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry diffDirAbs: ",ok) + print("BDTMuonEnergyPredict Tool: type of diffDirAbs entry is :",EnergyRecoBoostStore.Type("diffDirAbs")) + print("BDTMuonEnergyPredict Tool: Getting diffDirAbs from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("diffDirAbs",diffDirAbs) + print("BDTMuonEnergyPredict Tool: DiffDirAbs is: ", diffDirAbs.value) + ok = EnergyRecoBoostStore.Has("recoDWallR") + recoDWallR=ctypes.c_float(0) + if ok: + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry recoDWallR: ",ok) + print("BDTMuonEnergyPredict Tool: type of recoDWallR entry is :",EnergyRecoBoostStore.Type("recoDWallR")) + print("BDTMuonEnergyPredict Tool: Getting recoDWallR from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("recoDWallR",recoDWallR) + print("BDTMuonEnergyPredict Tool: RecoDWallR is: ", recoDWallR.value) + ok = EnergyRecoBoostStore.Has("recoDWallZ") + recoDWallZ=ctypes.c_float(0) + if(ok): + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry recoDWallZ: ",ok) + print("BDTMuonEnergyPredict Tool: type of recoDWallZ entry is :",EnergyRecoBoostStore.Type("recoDWallZ")) + print("BDTMuonEnergyPredict Tool: Getting recoDWallZ from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("recoDWallZ",recoDWallZ) + print("BDTMuonEnergyPredict Tool: RecoDWallZ is: ", recoDWallZ.value) + ok = EnergyRecoBoostStore.Has("vtxVec") + vtx_position=cppyy.gbl.Position() + if ok: + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry vtxVec: ",ok) + print("BDTMuonEnergyPredict Tool: type of vtxVec entry is :", EnergyRecoBoostStore.Type("vtxVec")) + print("BDTMuonEnergyPredict Tool: Getting vtxVec from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("vtxVec", vtx_position) + vtxX=vtx_position.X() + print("BDTMuonEnergyPredict Tool: VtxX is: ", vtxX) + vtxY=vtx_position.Y() + print("BDTMuonEnergyPredict Tool: VtxY is: ", vtxY) + vtxZ=vtx_position.Z() + print("BDTMuonEnergyPredict Tool: VtxZ is: ", vtxZ) + ok = EnergyRecoBoostStore.Has("trueMuonEnergy") + trueMuonEnergy=ctypes.c_double(0) + if ok: + print("BDTMuonEnergyPredict Tool: EnergyRecoBoostStore has entry trueMuonEnergy: ",ok) + print("BDTMuonEnergyPredict Tool: type of trueMuonEnergy entry is :",EnergyRecoBoostStore.Type("trueMuonEnergy")) + print("BDTMuonEnergyPredict Tool: Getting trueMuonEnergy from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("trueMuonEnergy",trueMuonEnergy) + print("BDTMuonEnergyPredict Tool: The MC muon energy is: ", trueMuonEnergy.value) + + #Create features and labels and preprocess data for the model + features_list=[] + features_list.append(DNNRecoLength.value/600.) + features_list.append(recoTrackLengthInMrd.value/200.) + features_list.append(diffDirAbs.value) + features_list.append(recoDWallR.value) + features_list.append(recoDWallZ.value) + features_list.append(num_lappd_hits.value/200.) + features_list.append(num_pmt_hits.value/200.) + features_list.append(vtxX/150.) + features_list.append(vtxY/200.) + features_list.append(vtxZ/150.) + features=np.array(features_list).reshape(1,10) + labels=np.array([trueMuonEnergy.value]) + + # load the model from disk + print("BDTMuonEnergyPredict Tool: Loading model") + loaded_model = pickle.load(open(str(self.weightsfilename), 'rb')) + + #predicting... + print("BDTMuonEnergyPredict Tool: Predicting") + recoEnergy = loaded_model.predict(features) + + #Set the BDTMuonEnergy in the EnergyReco boost store to be loaded by other tools + BDTMuonEnergy=float(recoEnergy[0]) + EnergyRecoBoostStore.Set("BDTMuonEnergy", BDTMuonEnergy) + self.m_data.Stores.at("ANNIEEvent").Set("RecoMuonEnergy", BDTMuonEnergy/1000.)#Set the BDTMuonEnergy in the ANNIEEvent store in [GeV] + + EnergyRecoBoostStore.Save("EnergyRecoStore.bs")#Append each entry to file in case we want to make plots + + return 1 + + def Finalise(self): + self.m_log.Log(__file__+" Finalising", self.v_debug, self.m_verbosity) + return 1 + +################### +# ↓ Boilerplate ↓ # +################### + +thistool = BDTMuonEnergyPredict() + +def SetToolChainVars(m_data_in, m_variables_in, m_log_in): + return thistool.SetToolChainVars(m_data_in, m_variables_in, m_log_in) + +def Initialise(): + return thistool.Initialise() + +def Execute(): + return thistool.Execute() + +def Finalise(): + return thistool.Finalise() diff --git a/UserTools/BDTMuonEnergyPredict/README.md b/UserTools/BDTMuonEnergyPredict/README.md new file mode 100644 index 000000000..390cbdba7 --- /dev/null +++ b/UserTools/BDTMuonEnergyPredict/README.md @@ -0,0 +1,40 @@ +# BDTMuonEnergyPredict + +The `BDTMuonEnergyPredict` tool is used to reconstruct the energy of the muon using a Boosted Decision Tree(BDT) with Gradient Boost. The tool loads the necessary information from the `EnergyReco` store. The `FindTrackLengthInWater` and `DNNTrackLengthPredict` tools should always be placed before the `BDTMuonEnergyPredict` tool in a toolchain. This tool loads the weights file of the model which has already been trained. The value of the reconstructed muon energy is passed to the next tools through the `EnergyReco` store. Also, this value is set into the `ANNIEEvent` store with the key `RecoMuonEnergy`. + +## Data + +### Input + +The following variables are obtained from the `EnergyReco` store: + +**DNNRecoLength** `double` The reconstructed track length in water for a single event + +**num_pmt_hits** `int` Total number of pmt digits + +**num_lappd_hits** `int` Total number of lappd digits + +**recoTrackLengthInMrd** Track length of a reconstructed track in the MRD found by the `FindMrdTracks` tool + +**diffDirAbs** `double` Angle difference between the reconstructed z direction and the beam direction at (0,0,1) + +**recoDWallR** `double` Radial distance of the reconstructed vertex from the walls of the tank + +**recoDWallZ** `double` Axial distance of the reconstructed vertex from the walls of the tank + +**vtxVec** `Position` Position of the reconstructed vertex + +**trueMuonEnergy** `double` MC muon energy + +### Output + +The following variables are passed on to the next tool via the `EnergyReco` store: + +**BDTMuonEnergy** `double` The reconstructed energy of the muon for a single event + +## Configuration + +``` +BDTMuonModelFile finalized_BDTmodel_forMuonEnergy.sav +The path of the weights file as it has been set when training the model +``` diff --git a/UserTools/BDTMuonEnergyTrain_Test/BDTMuonEnergyTrain_Test.py b/UserTools/BDTMuonEnergyTrain_Test/BDTMuonEnergyTrain_Test.py new file mode 100644 index 000000000..ad26fe60f --- /dev/null +++ b/UserTools/BDTMuonEnergyTrain_Test/BDTMuonEnergyTrain_Test.py @@ -0,0 +1,280 @@ +# BDTMuonEnergyTrain_Test Tool script# ------------------ +import sys +import numpy as np +import pandas as pd +import tempfile +import random +import csv +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from array import array +from sklearn import datasets +from sklearn import metrics +from sklearn import model_selection +from sklearn import preprocessing +import sklearn +from sklearn.utils import shuffle +from sklearn import linear_model, ensemble +from sklearn.metrics import mean_squared_error +import pickle +from Tool import * + +class BDTMuonEnergyTrain_Test(Tool): + + # declare member variables here + weightsfilename=std.string("") + inputBoostStorepathname = std.string("") + + def Initialise(self): + self.m_log.Log(__file__+" Initialising", self.v_debug, self.m_verbosity) + self.m_variables.Print() + self.m_variables.Get("BDTMuonEnergyWeightsFile", self.weightsfilename) + self.m_variables.Get("BDTMuonEnergyTrainingInputBoostStoreFile", self.inputBoostStorepathname) + return 1 + + def Execute(self): + self.m_log.Log(__file__+" Executing", self.v_debug, self.m_verbosity) + + #Set seed for reproducibility + seed=150 + + EnergyRecoBoostStore=cppyy.gbl.BoostStore(True, 2)#define the energy boost store class object to load the variables for the BDT training + ok=EnergyRecoBoostStore.Initialise(self.inputBoostStorepathname)#read from disk + print("BDTMuonEnergyTrain_Test Tool: Initiliased boost store successfully",ok) + total_entries = ctypes.c_ulong(0) + get_ok = EnergyRecoBoostStore.Header.Get("TotalEntries",total_entries) + print("BDTMuonEnergyTrain_Test Tool: Get num of entries of Energy Reco Store: ",get_ok,", entries: ",total_entries.value) + ievt=ctypes.c_ulong(0) + while True: + get_ok=EnergyRecoBoostStore.GetEntry(ievt.value) + print("BDTMuonEnergyTrain_Test Tool: There is an entry in the BoostStore",get_ok) + if not get_ok: + break; + #When there is no other entry GetEntry() returns false so the while loop stops + #Retrieve the required variables from this entry + EnergyRecoBoostStore.Print(False) + ok = EnergyRecoBoostStore.Has("num_pmt_hits") + num_pmt_hits=ctypes.c_int(0) + if ok: + print("BDTMuonEnergyTrain_Test Tool: EnergyRecoBoostStore has entry num_pmt_hits: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of num_pmt_hits entry is :",EnergyRecoBoostStore.Type("num_pmt_hits")) + print("BDTMuonEnergyTrain_Test Tool: Getting num_pmt_hits from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("num_pmt_hits",num_pmt_hits) + print("BDTMuonEnergyTrain_Test Tool: Number of pmt hits is: ", num_pmt_hits.value) + ok = EnergyRecoBoostStore.Has("num_lappd_hits") + num_lappd_hits=ctypes.c_int(0) + if ok: + print("BDTMuonEnergyTrain_Test Tool: EnergyRecoBoostStore has entry num_lappd_hits: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of num_lappd_hits entry is :",EnergyRecoBoostStore.Type("num_lappd_hits")) + print("BDTMuonEnergyTrain_Test Tool: Getting num_lappd_hits from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("num_lappd_hits",num_lappd_hits) + print("BDTMuonEnergyTrain_Test Tool: Number of lappd hits is: ", num_lappd_hits.value) + ok = EnergyRecoBoostStore.Has("recoTrackLengthInMrd") + recoTrackLengthInMrd=ctypes.c_double(0) + if ok: + print("BDTMuonEnergyTrain_Test Tool: EnergyRecoBoostStore has entry recoTrackLengthInMrd: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of recoTrackLengthInMrd entry is :",EnergyRecoBoostStore.Type("recoTrackLengthInMrd")) + print("BDTMuonEnergyTrain_Test Tool: Getting recoTrackLengthInMrd from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("recoTrackLengthInMrd",recoTrackLengthInMrd) + print("BDTMuonEnergyTrain_Test Tool: The reconstructed track length in the MRD is: ", recoTrackLengthInMrd.value) + ok = EnergyRecoBoostStore.Has("diffDirAbs") + diffDirAbs=ctypes.c_float(0) + if ok: + print("BDTMuonEnergyTrain_Test Tool: EnergyRecoBoostStore has entry diffDirAbs: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of diffDirAbs entry is :",EnergyRecoBoostStore.Type("diffDirAbs")) + print("BDTMuonEnergyTrain_Test Tool: Getting diffDirAbs from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("diffDirAbs",diffDirAbs) + print("BDTMuonEnergyTrain_Test Tool: DiffDirAbs is: ", diffDirAbs.value) + ok = EnergyRecoBoostStore.Has("recoDWallR") + recoDWallR=ctypes.c_float(0) + if ok: + print("BDTMuonEnergyTrain_Test Tool: EnergyRecoBoostStore has entry recoDWallR: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of recoDWallR entry is :",EnergyRecoBoostStore.Type("recoDWallR")) + print("BDTMuonEnergyTrain_Test Tool: Getting recoDWallR from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("recoDWallR",recoDWallR) + print("BDTMuonEnergyTrain_Test Tool: RecoDWallR is: ", recoDWallR.value) + ok = EnergyRecoBoostStore.Has("recoDWallZ") + recoDWallZ=ctypes.c_float(0) + if ok: + print("BDTMuonEnergyTrain_Test Tool: EnergyRecoBoostStore has entry recoDWallZ: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of recoDWallZ entry is :",EnergyRecoBoostStore.Type("recoDWallZ")) + print("BDTMuonEnergyTrain_Test Tool: Getting recoDWallZ from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("recoDWallZ",recoDWallZ) + print("BDTMuonEnergyTrain_Test Tool: RecoDWallZ is: ", recoDWallZ.value) + ok = EnergyRecoBoostStore.Has("vtxVec") + vtx_position=cppyy.gbl.Position() + if ok: + print("BDTMuonEnergyTrain_Test Tool: EnergyRecoBoostStore has entry vtxVec: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of vtxVec entry is :", EnergyRecoBoostStore.Type("vtxVec")) + print("BDTMuonEnergyTrain_Test Tool: Getting vtxVec from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("vtxVec", vtx_position) + vtxX=vtx_position.X() + print("BDTMuonEnergyTrain_Test Tool: VtxX is: ", vtxX) + vtxY=vtx_position.Y() + print("BDTMuonEnergyTrain_Test Tool: VtxY is: ", vtxY) + vtxZ=vtx_position.Z() + print("BDTMuonEnergyTrain_Test Tool: VtxZ is: ", vtxZ) + ok = EnergyRecoBoostStore.Has("trueMuonEnergy") + trueMuonEnergy=ctypes.c_double(0) + if ok: + print("BDTMuonEnergyTrain_Test Tool: EnergyRecoBoostStore has entry trueMuonEnergy: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of trueMuonEnergy entry is :",EnergyRecoBoostStore.Type("trueMuonEnergy")) + print("BDTMuonEnergyTrain_Test Tool: Getting trueMuonEnergy from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("trueMuonEnergy",trueMuonEnergy) + print("BDTMuonEnergyTrain_Test Tool: The MC muon energy is: ", trueMuonEnergy.value) + #Create features and labels and preprocess data for the model + features_list=[] + features_list.append(recoTrackLengthInMrd.value/200.) + features_list.append(diffDirAbs.value) + features_list.append(recoDWallR.value) + features_list.append(recoDWallZ.value) + features_list.append(num_lappd_hits.value/200.) + features_list.append(num_pmt_hits.value/200.) + features_list.append(vtxX/150.) + features_list.append(vtxY/200.) + features_list.append(vtxZ/150.) + #make the features and labels numpy array for this entry + featuresforthisentry=np.array(features_list) + labelsforthisentry=np.array([trueMuonEnergy.value]) + #vstack each entry + if ievt.value==0: + features=featuresforthisentry + labels=labelsforthisentry + else: + features=np.vstack([features,featuresforthisentry]) + labels=np.vstack([labels,labelsforthisentry]) + ievt.value+=1 + + print(features) + print(features.shape) + print(labels) + print(labels.shape) + + Dataset=np.concatenate((features,labels),axis=1) + + print(Dataset) + + np.random.seed(seed) + np.random.shuffle(Dataset) + + print(Dataset) + + features, labels = np.split(Dataset,[9],axis=1) + + num_events, num_pixels = features.shape + + print(num_events, num_pixels) + np.random.seed(0) + #select the events for which we have the reconstructed track length from the DNN + if len(labels) % 2==0: + a=int(len(labels)/2) + else: + a=int((len(labels)+1)/2) + features= features[a:] + labels = labels[a:] + + num_events, num_pixels = features.shape + print("BDTMuonEnergyTrain_Test Tool: Events with DNNRecoLength: ", num_events, num_pixels) + + #Get the DNN reconstructed track length in the water tank from RecoLength boost store + #Get the RecoLength Boost Store from DataModel + DNNRecoLengthBoostStore = self.m_data.Stores.at("RecoLength") + ok = DNNRecoLengthBoostStore.Has("DNNRecoLength") + DNNRecoLength_vector=std.vector[float](range(num_events)) + if ok: + print("BDTMuonEnergyTrain_Test Tool: DNNRecoLengthBoostStore has entry DNNRecoLength: ",ok) + print("BDTMuonEnergyTrain_Test Tool: type of DNNRecoLength entry is :", DNNRecoLengthBoostStore.Type("DNNRecoLength")) + print("BDTMuonEnergyTrain_Test Tool: Getting DNNRecoLength from DNNRecoLengthBoostStore") + DNNRecoLengthBoostStore.Get("DNNRecoLength", DNNRecoLength_vector) + print("BDTMuonEnergyTrain_Test Tool: The DNNRecoLength for the first event: ", DNNRecoLength_vector.at(0)) + DNNRecoLength_list=[] + for i in range(DNNRecoLength_vector.size()): + DNNRecoLength_list.append(DNNRecoLength_vector.at(i)/600.) + DNNRecoLength=np.array(DNNRecoLength_list).reshape(len(DNNRecoLength_list),1) + #Place the DNNRecoLength in the dataset + features=np.concatenate((DNNRecoLength, features), axis=1) + print(features) + print("BDTMuonEnergyTrain_Test Tool: Features shape after adding DNNRecoLength",features.shape) + + DNNRecoLength_new=features[:,0] + #This loop excludes any events with reconstructed length >1000 as not well reconstructed + j=0 + b=[] + for k in DNNRecoLength_new: + if k>1000: + print("RecoLength:",k,"Event:",j) + b.append(j) + j+=1 + features=np.delete(features,b,axis=0) + labels=np.delete(labels,b,axis=0) + + ########### BDTG ############ + n_estimators=500 + params = {'n_estimators':n_estimators, 'max_depth': 50, + 'learning_rate': 0.025, 'loss': 'absolute_error'} + + #--- select 70% of sample for training and 30% for testing: + offset = int(features.shape[0] * 0.7) + train_X = features[:offset] # train sample + train_Y = labels[:offset].reshape(-1)#reshape for model + test_X = features[offset:] # test sample + test_Y = labels[offset:].reshape(-1)#reshape for model + + print("BDTMuonEnergyTrain_Test Tool: train shape: ", train_X.shape," label: ", train_Y.shape) + print("BDTMuonEnergyTrain_Test Tool: test shape: ", test_X.shape," label: ", test_Y.shape) + + print("BDTMuonEnergyTrain_Test Tool: Training BDTG...") + net_hi_E = ensemble.GradientBoostingRegressor(**params) + model = net_hi_E.fit(train_X, train_Y) + net_hi_E + + # save the model to disk + filename = str(self.weightsfilename) + pickle.dump(model, open(filename, 'wb')) + + mse = mean_squared_error(test_Y, net_hi_E.predict(test_X)) + print("BDTMuonEnergyTrain_Test Tool: MSE: %.4f" % mse) + print("BDTMuonEnergyTrain_Test Tool: events at training & test samples: ", len(labels)) + print("BDTMuonEnergyTrain_Test Tool: events at train sample: ", len(train_Y)) + print("BDTMuonEnergyTrain_Test Tool: events at test sample: ", len(test_Y)) + + test_score = np.zeros((params['n_estimators'],), dtype=np.float64) + + for i, y_pred in enumerate(net_hi_E.staged_predict(test_X)): + test_score[i] = net_hi_E.loss_(test_Y, y_pred) + + fig,ax=plt.subplots(ncols=1, sharey=True) + ax.plot(np.arange(params['n_estimators']) + 1, net_hi_E.train_score_, 'b-', label='Training Set Error') + ax.plot(np.arange(params['n_estimators']) + 1, test_score, 'r-', label='Test Set Error') + ax.set_ylim(0.,500.) + ax.set_xlim(0.,n_estimators) + ax.legend(loc='upper right') + ax.set_ylabel('Absolute Errors [MeV]') + ax.set_xlabel('Number of Estimators') + ax.yaxis.set_label_coords(-0.1, 0.6) + ax.xaxis.set_label_coords(0.85, -0.08) + plt.savefig("error_train_test.png") + return 1 + + def Finalise(self): + self.m_log.Log(__file__+" Finalising", self.v_debug, self.m_verbosity) + return 1 + +################### +# ↓ Boilerplate ↓ # +################### + +thistool = BDTMuonEnergyTrain_Test() + +def SetToolChainVars(m_data_in, m_variables_in, m_log_in): + return thistool.SetToolChainVars(m_data_in, m_variables_in, m_log_in) + +def Initialise(): + return thistool.Initialise() + +def Execute(): + return thistool.Execute() + +def Finalise(): + return thistool.Finalise() diff --git a/UserTools/BDTMuonEnergyTrain_Test/README.md b/UserTools/BDTMuonEnergyTrain_Test/README.md new file mode 100644 index 000000000..c36050f66 --- /dev/null +++ b/UserTools/BDTMuonEnergyTrain_Test/README.md @@ -0,0 +1,44 @@ +# BDTMuonEnergyTrain_Test + +The `BDTMuonEnergyTrain_Test` tool is used to train a Boosted Decision Tree(BDT) with Gradient Boost which will then be able to reconstruct the energy of the muon. The tool loads the necessary information from the `EnergyReco` store which have been saved to disk by the `FindTrackLengthInWater` tool. The `FindTrackLengthInWater` tool should hence always be placed before the `BDTMuonEnergyTrain_Test` tool in a toolchain. The tool also loads the values for the reconstructed track length in the water for the each of the events which are calculated by the `DNNTrackLengthTrain_Test` tool and are passed to the `BDTMuonEnergyTrain_Test` tool through the `DNNRecoLength` store. The `DNNTrackLengthTrain_Test` tool should also always be placed before the `BDTMuonEnergyTrain_Test` tool in a toolchain. After training and testing the weights are saved locally. A few plots are also drawn for visualization purposes. + +## Data + +### Input + +The following variables are obtained from the `EnergyReco` store: + +**num_pmt_hits** `int` Total number of pmt digits + +**num_lappd_hits** `int` Total number of lappd digits + +**recoTrackLengthInMrd** Track length of a reconstructed track in the MRD found by the `FindMrdTracks` tool + +**diffDirAbs** `double` Angle difference between the reconstructed z direction and the beam direction at (0,0,1) + +**recoDWallR** `double` Radial distance of the reconstructed vertex from the walls of the tank + +**recoDWallZ** `double` Axial distance of the reconstructed vertex from the walls of the tank + +**vtxVec** `Position` Position of the reconstructed vertex + +**trueMuonEnergy** `double` MC muon energy + +The following variables are obtained from the `DNNRecoLength` store: + +**DNNRecoLength** `std::vector` Vector with the reconstructed track length in water for each event of the event sample + +## Configuration + +``` +InitialiseFunction Initialise +ExecuteFunction Finalise +FinaliseFunction Execute +Everything is done in the Execute method of the tool. We need to run the Execute method in the Finalise step of the toolchain so that the FindTrackLengthInWater tool has already saved a multiple event sample for the training. + +BDTMuonEnergyWeightsFile finalized_BDTmodel_forMuonEnergy.sav +The path where you want to save the weights file + +BDTMuonEnergyTrainingInputBoostStoreFile EnergyReco.bs +The path of the EnergyReco boost store +``` diff --git a/UserTools/DNNTrackLengthPredict/DNNTrackLengthPredict.py b/UserTools/DNNTrackLengthPredict/DNNTrackLengthPredict.py new file mode 100644 index 000000000..aef8b11ef --- /dev/null +++ b/UserTools/DNNTrackLengthPredict/DNNTrackLengthPredict.py @@ -0,0 +1,201 @@ +##### DNNTrackLengthPredict Tool Script +import numpy +import tensorflow +import random +import sys +import glob +import numpy as np +import pandas #as pd +import tempfile +import csv +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot #as plt +from array import array +from sklearn import datasets +from sklearn import metrics +from sklearn import model_selection +from sklearn import preprocessing +from tensorflow import keras +from tensorflow.keras.models import Sequential +from tensorflow.keras.layers import Dense +from tensorflow.keras.callbacks import ModelCheckpoint +from tensorflow.keras.wrappers.scikit_learn import KerasRegressor +from tensorflow.keras import backend as K +import pprint +#import Store and other modules +from Tool import * + +class DNNTrackLengthPredict(Tool): + + # declare member variables here + weightsfilename = std.string("") + ScalingVarsBoostStorepathname = std.string("") + + def Initialise(self): + self.m_log.Log(__file__+" Initialising", self.v_debug, self.m_verbosity) + self.m_variables.Print() + self.m_variables.Get("TrackLengthWeightsFile", self.weightsfilename) + self.m_variables.Get("ScalingVarsBoostStoreFile", self.ScalingVarsBoostStorepathname) + return 1 + + def Execute(self): + self.m_log.Log(__file__+" Executing", self.v_debug, self.m_verbosity) + # Load Data from EnergyReco Boost Store + #----------------------------- + EnergyRecoBoostStore=cppyy.gbl.BoostStore(True, 2)#define the energy boost store class object to load the variables for the DNN training + EnergyRecoBoostStore = self.m_data.Stores.at("EnergyReco") + #Retrieve the required variables from the Store + EnergyRecoBoostStore.Print(False) + get_ok = EnergyRecoBoostStore.Has("MaxTotalHitsToDNN") + MaxTotalHitsToDNN=ctypes.c_int(0) + if not get_ok: + print("DNNTrackLengthPredict Tool: There is no entry in Energy Reco boost store.") + return 1 + if get_ok: + print("DNNTrackLengthPredict Tool: EnergyRecoBoostStore has entry MaxTotalHitsToDNN: ",get_ok) + print("DNNTrackLengthPredict Tool: type of MaxTotalHitsToDNN entry is :",EnergyRecoBoostStore.Type("MaxTotalHitsToDNN")) + print("DNNTrackLengthPredict Tool: Getting MaxTotalHitsToDNN from EnergyRecoBoostStore")#we are going to use it to instantiate the lambda and digit times vectors + EnergyRecoBoostStore.Get("MaxTotalHitsToDNN",MaxTotalHitsToDNN) + print("DNNTrackLengthPredict Tool: MaxTotalHitsToDNN is: ", MaxTotalHitsToDNN.value) + ok = EnergyRecoBoostStore.Has("lambda_vec") + lambda_vector=std.vector['double'](range(MaxTotalHitsToDNN.value)) + if ok: + print("DNNTrackLengthPredict Tool: EnergyRecoBoostStore has entry lambda_vec: ",ok) + print("DNNTrackLengthPredict Tool: type of lambda_vec entry is :", EnergyRecoBoostStore.Type("lambda_vec")) + print("DNNTrackLengthPredict Tool: Getting lambda_vec from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("lambda_vec", lambda_vector) + print("DNNTrackLengthPredict Tool: The lambda for the first digit is: ", lambda_vector.at(0)) + ok = EnergyRecoBoostStore.Has("digit_ts_vec") + digit_ts_vector=std.vector['double'](range(MaxTotalHitsToDNN.value)) + if ok: + print("DNNTrackLengthPredict Tool: EnergyRecoBoostStore has entry digit_ts_vec: ",ok) + print("DNNTrackLengthPredict Tool: type of digit_ts_vec entry is :",EnergyRecoBoostStore.Type("digit_ts_vec")) + print("DNNTrackLengthPredict Tool: Getting digit_ts_vec from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("digit_ts_vec", digit_ts_vector) + print("DNNTrackLengthPredict Tool: The digit time for the first digit is: ", digit_ts_vector.at(0)) + ok = EnergyRecoBoostStore.Has("lambda_max") + lambda_max=ctypes.c_double(0) + if ok: + print("DNNTrackLengthPredict Tool: EnergyRecoBoostStore has entry lambda_max: ",ok) + print("DNNTrackLengthPredict Tool: type of lambda_max entry is :",EnergyRecoBoostStore.Type("lambda_max")) + print("DNNTrackLengthPredict Tool: Getting lambda_max from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("lambda_max",lambda_max) + print("DNNTrackLengthPredict Tool: Lambda_max is: ", lambda_max.value) + ok = EnergyRecoBoostStore.Has("num_pmt_hits") + num_pmt_hits=ctypes.c_int(0) + if ok: + print("DNNTrackLengthPredict Tool: EnergyRecoBoostStore has entry num_pmt_hits: ",ok) + print("DNNTrackLengthPredict Tool: type of num_pmt_hits entry is :",EnergyRecoBoostStore.Type("num_pmt_hits")) + print("DNNTrackLengthPredict Tool: Getting num_pmt_hits from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("num_pmt_hits",num_pmt_hits) + print("DNNTrackLengthPredict Tool: Number of pmt hits is: ", num_pmt_hits.value) + ok = EnergyRecoBoostStore.Has("num_lappd_hits") + num_lappd_hits=ctypes.c_int(0) + if ok: + print("DNNTrackLengthPredict Tool: EnergyRecoBoostStore has entry num_lappd_hits: ",ok) + print("DNNTrackLengthPredict Tool: type of num_lappd_hits entry is :",EnergyRecoBoostStore.Type("num_lappd_hits")) + print("DNNTrackLengthPredict Tool: Getting num_lappd_hits from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("num_lappd_hits",num_lappd_hits) + print("DNNTrackLengthPredict Tool: Number of lappd hits is: ", num_lappd_hits.value) + ok = EnergyRecoBoostStore.Has("TrueTrackLengthInWater") + TrueTrackLengthInWater=ctypes.c_float(0) + if ok: + print("DNNTrackLengthPredict Tool: EnergyRecoBoostStore has entry TrueTrackLengthInWater: ",ok) + print("DNNTrackLengthPredict Tool: type of TrueTrackLengthInWater entry is :",EnergyRecoBoostStore.Type("TrueTrackLengthInWater")) + print("DNNTrackLengthPredict Tool: Getting TrueTrackLengthInWater from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("TrueTrackLengthInWater",TrueTrackLengthInWater) + print("DNNTrackLengthPredict Tool: TrueTrackLengthInWater is: ", TrueTrackLengthInWater.value) + + #Create features and labels and preprocess data for the model + features_list=[] + for i in range(lambda_vector.size()): + features_list.append(lambda_vector.at(i)) + for j in range(digit_ts_vector.size()): + features_list.append(digit_ts_vector.at(j)) + features_list.append(lambda_max.value) + features_list.append(num_pmt_hits.value) + features_list.append(num_lappd_hits.value) + #make the features and labels numpy array + features=np.array(features_list) + labels=np.array([TrueTrackLengthInWater.value]) + + print(features) + print(features.shape) + print(labels) + print(labels.shape) + + #Load scaling parameters + ScalingVarsStore=cppyy.gbl.BoostStore(True, 0) + ok=ScalingVarsStore.Initialise(self.ScalingVarsBoostStorepathname) + features_mean_values_vec=std.vector['double'](range(len(features))) + features_std_values_vec=std.vector['double'](range(len(features))) + ScalingVarsStore.Get("features_mean_values",features_mean_values_vec) + ScalingVarsStore.Get("features_std_values",features_std_values_vec) + + #Transform data + test_x=[] + + for i in range(len(features)): + test_x.append((features[i]-features_mean_values_vec.at(i))/features_std_values_vec.at(i)) + test_X=np.array(test_x).reshape(1,2203) + + print("DNNTrackLengthPredict Tool: Defining the model") + model = Sequential() + print("DNNTrackLengthPredict Tool: Adding layers") + model.add(Dense(25, input_dim=2203, kernel_initializer='normal', activation='relu')) + model.add(Dense(25, kernel_initializer='normal', activation='relu')) + model.add(Dense(1, kernel_initializer='normal', activation='relu')) + + # load weights + print("DNNTrackLengthPredict Tool: Loading weights from file ", self.weightsfilename) + model.load_weights(str(self.weightsfilename)) + + # Compile model + print("DNNTrackLengthPredict Tool: Compiling model") + model.compile(loss='mean_squared_error', optimizer='Adamax', metrics=['accuracy']) + print("DNNTrackLengthPredict Tool: Created model and loaded weights from file", self.weightsfilename) + + # Score accuracy / Make predictions + #---------------------------------- + print('DNNTrackLengthPredict Tool: Predicting...') + y_predicted = model.predict(test_X) + print(y_predicted.shape) + + # estimate accuracy on dataset using loaded weights + print("DNNTrackLengthPredict Tool: Evalulating model on test") + scores = model.evaluate(test_X, labels, verbose=0) + print("%s: %.2f%%" % (model.metrics_names[1], scores[1]*100)) + + # Score with sklearn. + print("scoring sk mse") + score_sklearn = metrics.mean_squared_error(y_predicted, labels) + print('MSE (sklearn): {0:f}'.format(score_sklearn)) + + #Set the DNNRecoLength in the EnergyReco boost store for next tools + DNNRecoLength=float(y_predicted[0]) + EnergyRecoBoostStore.Set("DNNRecoLength", DNNRecoLength) + self.m_data.Stores.at("ANNIEEvent").Set("RecoTrackLengthInTank", DNNRecoLength/100.)#Set the DNNRecoLength in the ANNIEEvent store in [m] + return 1 + + def Finalise(self): + self.m_log.Log(__file__+" Finalising", self.v_debug, self.m_verbosity) + return 1 + +################### +# ↓ Boilerplate ↓ # +################### + +thistool = DNNTrackLengthPredict() + +def SetToolChainVars(m_data_in, m_variables_in, m_log_in): + return thistool.SetToolChainVars(m_data_in, m_variables_in, m_log_in) + +def Initialise(): + return thistool.Initialise() + +def Execute(): + return thistool.Execute() + +def Finalise(): + return thistool.Finalise() diff --git a/UserTools/DNNTrackLengthPredict/README.md b/UserTools/DNNTrackLengthPredict/README.md new file mode 100644 index 000000000..f77bddcff --- /dev/null +++ b/UserTools/DNNTrackLengthPredict/README.md @@ -0,0 +1,45 @@ +# DNNTrackLengthPredict + +The `DNNTrackLengthPredict` tool is used to reconstruct the track length of the muon in the water tank using a Deep Learning Neural Network(DNN) from the Keras library on top of Tensorflow. The tool loads the necessary information from the `EnergyReco` store. The `FindTrackLengthInWater` tool should always be placed before the `DNNTrackLengthPredict` tool in a toolchain. The scaling parameters are loaded from the `ScalingVarsStore` in order to preprocess our data for the model. This tool also loads the weights file of the model which has already been trained. The value of the reconstructed track length in the water tank is passed to the next tools through the `EnergyReco` store. Also, this value is set into the `ANNIEEvent` store with the key `RecoTrackLengthInTank`. + +## Data + +### Input + +The following variables are obtained from the `EnergyReco` store: + +**MaxTotalHitsToDNN** `int` + +**lambda_vec** `std::vector` Vector with all the lambda values for each event + +**digit_ts_vec** `std::vector` Vector with the time of all the digits of each event + +**lambda_max** `double` The distance between the reconstructed vertex and last Cherenkov photon emission point along the track + +**num_pmt_hits** `int` Total number of pmt digits + +**num_lappd_hits** `int` Total number of lappd digits + +**TrueTrackLengthInWater** `float` MC track length in the water + +The following variables are obtained from the `ScalingVarsStore` store: + +**features_mean_values** `std::vector` Vector with the means of each feature of the training dataset + +**features_std_values** `std::vector` Vector with the standard deviations of each feature of the training dataset + +### Output + +The following variables are passed on to the next tool via the `EnergyReco` store: + +**DNNRecoLength** `double` The reconstructed track length in water for a single event + +## Configuration + +``` +TrackLengthWeightsFile UserTools/DNNTrackLength/stand_alone/weights/weights_bets.hdf5 +The path of the weights file as it has been set when training the model + +ScalingVarsBoostStoreFile Data_Energy_Reco/ScalingVarsStore.bs +The path of the scaling parameters file as it has been set when training the model +``` diff --git a/UserTools/DNNTrackLengthTrain_Test/DNNTrackLengthTrain_Test.py b/UserTools/DNNTrackLengthTrain_Test/DNNTrackLengthTrain_Test.py new file mode 100644 index 000000000..65a678344 --- /dev/null +++ b/UserTools/DNNTrackLengthTrain_Test/DNNTrackLengthTrain_Test.py @@ -0,0 +1,335 @@ +# DNNTrackLengthTrain_Test Tool script +# ------------------ +import sys +import glob +import numpy as np +import pandas as pd +import tensorflow as tf +import tempfile +import random +import csv +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from array import array +from sklearn import datasets +from sklearn import metrics +from sklearn import model_selection +from sklearn import preprocessing +from sklearn.utils import shuffle +from tensorflow import keras +from tensorflow.keras.models import Sequential +from tensorflow.keras.layers import Dense +from tensorflow.keras.callbacks import ModelCheckpoint +from tensorflow.keras.wrappers.scikit_learn import KerasRegressor +import ROOT +import matplotlib.pylab as pylab +params = {'legend.fontsize': 'x-large', + 'figure.figsize': (9, 7), + 'axes.labelsize': 'x-large', + 'axes.titlesize':'x-large', + 'xtick.labelsize':'x-large', + 'ytick.labelsize':'x-large'} +pylab.rcParams.update(params) +from Tool import * + +class DNNTrackLengthTrain_Test(Tool): + + # declare member variables here + weightsfilename = std.string("") + inputBoostStorepathname = std.string("") + ScalingVarsBoostStorepathname = std.string("") + + def Initialise(self): + self.m_log.Log(__file__+" Initialising", self.v_debug, self.m_verbosity) + self.m_variables.Print() + self.m_variables.Get("TrackLengthOutputWeightsFile", self.weightsfilename) + self.m_variables.Get("TrackLengthTrainingInputBoostStoreFile", self.inputBoostStorepathname) + self.m_variables.Get("ScalingVarsBoostStoreFile", self.ScalingVarsBoostStorepathname) + return 1 + + def Execute(self): + self.m_log.Log(__file__+" Executing", self.v_debug, self.m_verbosity) + #Set seed for reproducibility + + seed=150 + + EnergyRecoBoostStore=cppyy.gbl.BoostStore(True, 2)#define the energy boost store class object to load the variables for the DNN training + ok=EnergyRecoBoostStore.Initialise(self.inputBoostStorepathname)#read from disk + print("DNNTrackLengthTrain_Test Tool: Initiliased boost store successfully") + total_entries = ctypes.c_ulong(0) + get_ok = EnergyRecoBoostStore.Header.Get("TotalEntries",total_entries) + print("DNNTrackLengthTrain_Test Tool: Get num of entries of Energy Reco Store: ",get_ok,", entries: ",total_entries.value) + ievt=ctypes.c_ulong(0) + while True: + get_ok=EnergyRecoBoostStore.GetEntry(ievt.value) + print("DNNTrackLengthTrain_Test Tool: There is an entry in the BoostStore",get_ok) + if not get_ok: + break; + #When there is no other entry GetEntry() returns false so the while loop stops + #Retrieve the required variables from this entry + EnergyRecoBoostStore.Print(False) + ok = EnergyRecoBoostStore.Has("MaxTotalHitsToDNN") + MaxTotalHitsToDNN=ctypes.c_int(0) + if ok: + print("DNNTrackLengthTrain_Test Tool: EnergyRecoBoostStore has entry MaxTotalHitsToDNN: ",ok) + print("DNNTrackLengthTrain_Test Tool: type of MaxTotalHitsToDNN entry is :",EnergyRecoBoostStore.Type("MaxTotalHitsToDNN")) + print("DNNTrackLengthTrain_Test Tool: Getting MaxTotalHitsToDNN from EnergyRecoBoostStore")#we are going to use it to instantiate the lambda and digit times vectors + EnergyRecoBoostStore.Get("MaxTotalHitsToDNN",MaxTotalHitsToDNN) + print("DNNTrackLengthTrain_Test Tool: MaxTotalHitsToDNN is: ", MaxTotalHitsToDNN.value) + ok = EnergyRecoBoostStore.Has("lambda_vec") + lambda_vector=std.vector['double'](range(MaxTotalHitsToDNN.value)) + if ok: + print("DNNTrackLengthTrain_Test Tool: EnergyRecoBoostStore has entry lambda_vec: ",ok) + print("DNNTrackLengthTrain_Test Tool: type of lambda_vec entry is :", EnergyRecoBoostStore.Type("lambda_vec")) + print("DNNTrackLengthTrain_Test Tool: Getting lambda_vec from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("lambda_vec", lambda_vector) + print("DNNTrackLengthTrain_Test Tool: The lambda for the first digit is: ", lambda_vector.at(0)) + ok = EnergyRecoBoostStore.Has("digit_ts_vec") + digit_ts_vector=std.vector['double'](range(MaxTotalHitsToDNN.value)) + if ok: + print("DNNTrackLengthTrain_Test Tool: EnergyRecoBoostStore has entry digit_ts_vec: ",ok) + print("DNNTrackLengthTrain_Test Tool: type of digit_ts_vec entry is :",EnergyRecoBoostStore.Type("digit_ts_vec")) + print("DNNTrackLengthTrain_Test Tool: Getting digit_ts_vec from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("digit_ts_vec", digit_ts_vector) + print("DNNTrackLengthTrain_Test Tool: The digit time for the first digit is: ", digit_ts_vector.at(0)) + ok = EnergyRecoBoostStore.Has("lambda_max") + lambda_max=ctypes.c_double(0) + if ok: + print("DNNTrackLengthTrain_Test Tool: EnergyRecoBoostStore has entry lambda_max: ",ok) + print("DNNTrackLengthTrain_Test Tool: type of lambda_max entry is :",EnergyRecoBoostStore.Type("lambda_max")) + print("DNNTrackLengthTrain_Test Tool: Getting lambda_max from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("lambda_max",lambda_max) + print("DNNTrackLengthTrain_Test Tool: Lambda_max is: ", lambda_max.value) + ok = EnergyRecoBoostStore.Has("num_pmt_hits") + num_pmt_hits=ctypes.c_int(0) + if ok: + print("DNNTrackLengthTrain_Test Tool: EnergyRecoBoostStore has entry num_pmt_hits: ",ok) + print("DNNTrackLengthTrain_Test Tool: type of num_pmt_hits entry is :",EnergyRecoBoostStore.Type("num_pmt_hits")) + print("Getting num_pmt_hits from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("num_pmt_hits",num_pmt_hits) + print("DNNTrackLengthTrain_Test Tool: Number of pmt hits is: ", num_pmt_hits.value) + ok = EnergyRecoBoostStore.Has("num_lappd_hits") + num_lappd_hits=ctypes.c_int(0) + if ok: + print("DNNTrackLengthTrain_Test Tool: EnergyRecoBoostStore has entry num_lappd_hits: ",ok) + print("DNNTrackLengthTrain_Test Tool: type of num_lappd_hits entry is :",EnergyRecoBoostStore.Type("num_lappd_hits")) + print("DNNTrackLengthTrain_Test Tool: Getting num_lappd_hits from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("num_lappd_hits",num_lappd_hits) + print("DNNTrackLengthTrain_Test Tool: Number of lappd hits is: ", num_lappd_hits.value) + ok = EnergyRecoBoostStore.Has("TrueTrackLengthInWater") + TrueTrackLengthInWater=ctypes.c_float(0) + if ok: + print("DNNTrackLengthTrain_Test Tool: EnergyRecoBoostStore has entry TrueTrackLengthInWater: ",ok) + print("DNNTrackLengthTrain_Test Tool: type of TrueTrackLengthInWater entry is :",EnergyRecoBoostStore.Type("TrueTrackLengthInWater")) + print("DNNTrackLengthTrain_Test Tool: Getting TrueTrackLengthInWater from EnergyRecoBoostStore") + EnergyRecoBoostStore.Get("TrueTrackLengthInWater",TrueTrackLengthInWater) + print("DNNTrackLengthTrain_Test Tool: TrueTrackLengthInWater is: ", TrueTrackLengthInWater.value) + + #Create features and labels and preprocess data for the model + features_list=[] + for i in range(lambda_vector.size()): + features_list.append(lambda_vector.at(i)) + for j in range(digit_ts_vector.size()): + features_list.append(digit_ts_vector.at(j)) + features_list.append(lambda_max.value) + features_list.append(num_pmt_hits.value) + features_list.append(num_lappd_hits.value) + #make the features and labels numpy array for this entry + featuresforthisentry=np.array(features_list) + labelsforthisentry=np.array([TrueTrackLengthInWater.value]) + #vstack each entry + if ievt.value==0: + features=featuresforthisentry + labels=labelsforthisentry + else: + features=np.vstack([features,featuresforthisentry]) + labels=np.vstack([labels,labelsforthisentry]) + ievt.value+=1 + + print(features) + print(features.shape) + print(labels) + print(labels.shape) + + Dataset=np.concatenate((features,labels),axis=1) + + print(Dataset) + + #shuffle the data in order to avoid any bias in training + np.random.seed(seed) + np.random.shuffle(Dataset) + + print(Dataset) + + features, labels = np.split(Dataset,[2203],axis=1) + + #split events in train/test samples: + + num_events, num_pixels = features.shape + print(num_events, num_pixels) + np.random.seed(0) + #split in half if the number of events is even or take one more event for training when odd + if len(labels) % 2==0: + a=int(len(labels)/2) + else: + a=int((len(labels)+1)/2) + train_x = features[:a] + train_y = labels[:a] + test_x = features[a:] + test_y = labels[a:] + + print("DNNTrackLengthTrain_Test Tool: train sample features shape: ", train_x.shape," DNNTrackLengthTrain_Test Tool: train sample label shape: ", train_y.shape) + print("DNNTrackLengthTrain_Test Tool: test sample features shape: ", test_x.shape," DNNTrackLengthTrain_Test Tool: test sample label shape: ", test_y.shape) + + # Scale data (training set) to 0 mean and unit standard deviation. + scaler = preprocessing.StandardScaler() + train_x = scaler.fit_transform(train_x) + + #Save the scaling parameters from fit_transform for later use in the predict toolchain + features_mean_values=scaler.mean_ + features_std_values=scaler.scale_ + + #Turn into std.vectors() to store in the ScalingVars BoostStore + features_mean_values_vec=std.vector['double'](range(len(features_mean_values))) + for i in range(features_mean_values_vec.size()): + features_mean_values_vec[i]=features_mean_values[i] + features_std_values_vec=std.vector['double'](range(len(features_std_values))) + for j in range(features_std_values_vec.size()): + features_std_values_vec[j]=features_std_values[j] + #Create a boost store to store these values + ScalingVarsStore=cppyy.gbl.BoostStore(True, 0) + print("DNNTrackLengthTrain_Test Tool: Constructed: ",type(ScalingVarsStore)) + ScalingVarsStore.Set("features_mean_values",features_mean_values_vec) + ScalingVarsStore.Set("features_std_values",features_std_values_vec) + #Save BoostStore locally for later use in EnergyRecoPredict toolchain + print("DNNTrackLengthTrain_Test Tool: Saving boost store with scaling parameters locally for later use") + ScalingVarsStore.Save(self.ScalingVarsBoostStorepathname) + def create_model(): + # create model + model = Sequential() + model.add(Dense(25, input_dim=2203, kernel_initializer='normal', activation='relu')) + model.add(Dense(25, kernel_initializer='normal', activation='relu')) + model.add(Dense(1, kernel_initializer='normal', activation='relu')) + # Compile model + model.compile(loss='mean_squared_error', optimizer='Adamax', metrics=['accuracy']) + return model + + estimator = KerasRegressor(build_fn=create_model, epochs=10, batch_size=2, verbose=0) + + # checkpoint + filepath=str(self.weightsfilename) + checkpoint = ModelCheckpoint(filepath, monitor='val_loss', verbose=1, save_best_only=True, save_weights_only=True, mode='auto') + callbacks_list = [checkpoint] + # Fit the model + + history = estimator.fit(train_x, train_y, validation_split=0.33, epochs=12, batch_size=1, callbacks=callbacks_list, verbose=0) + #----------------------------- + # summarize history for loss + fig, ax = plt.subplots(1,1) + ax.plot(history.history['loss']) + ax.plot(history.history['val_loss']) + ax.set_title('Model Loss') + ax.set_ylabel('Performance') + ax.set_xlabel('Epochs') + ax.set_xlim(0.,12.) + ax.legend(['training loss', 'validation loss'], loc='upper left') + plt.savefig("keras_train_test.pdf") + plt.close(fig) + + # create model to test the weights + model = Sequential() + model.add(Dense(25, input_dim=2203, kernel_initializer='normal', activation='relu')) + model.add(Dense(25, kernel_initializer='normal', activation='relu')) + model.add(Dense(1, kernel_initializer='normal', activation='relu')) + + # load weights + model.load_weights(filepath) + + # Compile model + model.compile(loss='mean_squared_error', optimizer='Adamax', metrics=['accuracy']) + print("DNNTrackLengthTrain_Test Tool: Created model and loaded weights from"+filepath+"for testing") + + ## Predict. + print('DNNTrackLengthTrain_Test Tool: Predicting...') + x_transformed = scaler.transform(test_x) + y_predicted = model.predict(x_transformed) + + scores = model.evaluate(x_transformed, test_y, verbose=0) + print("%s: %.2f%%" % (model.metrics_names[1], scores[1]*100)) + + # Score with sklearn. + score_sklearn = metrics.mean_squared_error(y_predicted, test_y) + print('MSE (sklearn): {0:f}'.format(score_sklearn)) + + #----------------------------- + print("shapes: ", test_y.shape, ", ", y_predicted.shape) + #Creating a BoostStore to store the DNN reconstructed track length in the water in the DataModel to be loaded for the BDT trainining + RecoLengthStore = cppyy.gbl.BoostStore(True,0) + DNNRecoLength=std.vector[float](range(len(y_predicted))) + for k in range(len(y_predicted)): + DNNRecoLength[k]=float(y_predicted[k]) + print(DNNRecoLength.at(0),y_predicted[0]) + RecoLengthStore.Set("DNNRecoLength", DNNRecoLength) + self.m_data.Stores["RecoLength"] = RecoLengthStore + + #Make some plots for visualization purposes + fig0, ax0 = plt.subplots() + ax0.scatter(test_y,y_predicted) + ax0.plot([test_y.min(),test_y.max()],[test_y.min(),test_y.max()],'k--',lw=3) + ax0.set_xlabel('MC Tracklength in water(cm)') + ax0.set_ylabel('Predicted Reconstucted Length from DNN(cm)') + fig0.savefig('recolength.png') + plt.close(fig0) + + data = abs(y_predicted-test_y) + lambda_max=test_x[:,[2200]].tolist() + dataprev = abs(lambda_max-test_y) + nbins=np.arange(0.,400.,5) + fig1,ax1=plt.subplots(ncols=1, sharey=False) + f0=ax1.hist(data, nbins, histtype='step', fill=False, color='blue',alpha=0.75) + f1=ax1.hist(dataprev, nbins, histtype='step', fill=False, color='red',alpha=0.75) + ax1.set_xlabel('$\Delta R = |L_{Reco}-L_{MC}|$ [cm]') + ax1.legend(('NEW','Previous')) + ax1.xaxis.set_ticks(np.arange(0., 425., 25)) + ax1.tick_params(axis='x', which='minor', bottom=False) + xmin, xmax = plt.xlim() + title = "mean = %.2f, std = %.2f, Prev: mean = %.2f, std = %.2f " % (data.mean(), data.std(),dataprev.mean(), dataprev.std()) + plt.title(title) + fig1.savefig('resol_distr2l_WCSim.png') + plt.close(fig1) + + canvas = ROOT.TCanvas() + canvas.cd(1) + th2f = ROOT.TH2F("True_RecoLength", "; MC Track Length [cm]; Reconstructed Track Length [cm]", 50, 0, 400., 50, 0., 400.) + for i in range(len(test_y)): + th2f.Fill(test_y[i], y_predicted[i]) + line = ROOT.TLine(0.,0.,400.,400.) + th2f.SetStats(0) + th2f.Draw("ColZ") + line.SetLineColor(2) + line.Draw("same") + canvas.SaveAs("MClength_newrecolength.png") + return 1 + + def Finalise(self): + self.m_log.Log(__file__+" Finalising", self.v_debug, self.m_verbosity) + return 1 + +################### +# ↓ Boilerplate ↓ # +################### + +thistool = DNNTrackLengthTrain_Test() + +def SetToolChainVars(m_data_in, m_variables_in, m_log_in): + return thistool.SetToolChainVars(m_data_in, m_variables_in, m_log_in) + +def Initialise(): + return thistool.Initialise() + +def Execute(): + return thistool.Execute() + +def Finalise(): + return thistool.Finalise() diff --git a/UserTools/DNNTrackLengthTrain_Test/README.md b/UserTools/DNNTrackLengthTrain_Test/README.md new file mode 100644 index 000000000..906f07cb5 --- /dev/null +++ b/UserTools/DNNTrackLengthTrain_Test/README.md @@ -0,0 +1,49 @@ +# DNNTrackLengthTrain_Test + +The `DNNTrackLengthTrain_Test` tool is used to train a Deep Learning Neural Network(DNN) using the Keras library on top of Tensorflow to reconstruct the track length of the muon in the water. The tool loads the necessary information from the `EnergyReco` store which have saved to disk by the `FindTrackLengthInWater` tool. The `FindTrackLengthInWater` tool should hence always be placed before the `DNNTrackLengthTrain_Test` tool in a toolchain. The input dataset is split into half in order to use one half for training and testing the model and the second half for reconstructing the track length in the water and use these events to train the BDT for the muon energy in the next tool. After training and testing the weights are saved locally. The values of the reconstructed track length in the water for the second half of the dataset are saved in the `DNNRecoLength` store to be loaded by the next tool. Also, some parameters concerning the training dataset that will need to be used in the prediction toolchain to scale our data are saved locally in the `ScalingVarsStore`. A few plots are also drawn for visualization purposes. + +## Data + +### Input + +The following variables are obtained from the `EnergyReco` store: + +**MaxTotalHitsToDNN** `int` + +**lambda_vec** `std::vector` Vector with all the lambda values for each event + +**digit_ts_vec** `std::vector` Vector with the time of all the digits of each event + +**lambda_max** `double` The distance between the reconstructed vertex and last Cherenkov photon emission point along the track + +**num_pmt_hits** `int` Total number of pmt digits + +**num_lappd_hits** `int` Total number of lappd digits + +**TrueTrackLengthInWater** `float` MC track length in the water + +### Output + +**features_mean_values** `std::vector` Vector with the means of each feature of the training dataset + +**features_std_values** `std::vector` Vector with the standard deviations of each feature of the training dataset + +**DNNRecoLength** `std::vector` Vector with the reconstructed track length in water for each event of the second half of the event sample + +## Configuration + +``` +InitialiseFunction Initialise +ExecuteFunction Finalise +FinaliseFunction Execute +Everything is done in the Execute method of the tool. We need to run the Execute method in the Finalise step of the toolchain so that the FindTrackLengthInWater tool has already saved a multiple event sample for the training. + +TrackLengthOutputWeightsFile UserTools/DNNTrackLength/stand_alone/weights/weights_bets.hdf5 +The path where you want to save the weights file + +TrackLengthTrainingInputBoostStoreFile EnergyReco.bs +The path of the EnergyReco boost store + +ScalingVarsBoostStoreFile Data_Energy_Reco/ScalingVarsStore.bs +The path where you want to save the scaling parameters file +``` diff --git a/UserTools/Factory/Factory.cpp b/UserTools/Factory/Factory.cpp index 2215596e1..cd6659bd5 100644 --- a/UserTools/Factory/Factory.cpp +++ b/UserTools/Factory/Factory.cpp @@ -159,6 +159,7 @@ if (tool=="FilterLAPPDEvents") ret=new FilterLAPPDEvents; if (tool=="FilterEvents") ret=new FilterEvents; if (tool=="Stage1DataBuilder") ret=new Stage1DataBuilder; if (tool=="BeamFetcherV2") ret=new BeamFetcherV2; +if (tool=="PlotsTrackLengthAndEnergy") ret=new PlotsTrackLengthAndEnergy; if (tool=="SaveConfigInfo") ret=new SaveConfigInfo; if (tool=="ReadConfigInfo") ret=new ReadConfigInfo; return ret; diff --git a/UserTools/FindTrackLengthInWater/FindTrackLengthInWater.cpp b/UserTools/FindTrackLengthInWater/FindTrackLengthInWater.cpp index 260fde121..205724970 100644 --- a/UserTools/FindTrackLengthInWater/FindTrackLengthInWater.cpp +++ b/UserTools/FindTrackLengthInWater/FindTrackLengthInWater.cpp @@ -1,4 +1,6 @@ #include "FindTrackLengthInWater.h" +#include +#include "TMath.h" FindTrackLengthInWater::FindTrackLengthInWater():Tool(){} @@ -11,226 +13,357 @@ bool FindTrackLengthInWater::Initialise(std::string configfile, DataModel &data) m_data= &data; //assigning transient data pointer ///////////////////////////////////////////////////////////////// - // get configuration variables for this tool - m_variables.Get("InputFile",infile); - - file= new TFile(infile.c_str(),"READ"); - regTree= (TTree*) file->Get("vertextree"); - std::cout<<"Number of entries in tree: "<GetEntries()<1100){ - std::cerr<<" Please change the dim of double lambda_vec[1100]={0.}; double digitt[1100]={0.}; from 1100 to max number of hits"<Stores.emplace("EnergyReco",energystore);// place boost store in the DataModel + + // Get values from Config file + // =========================== + get_ok = m_variables.Get("MaxTotalHitsToDNN",maxhits0); + if(not get_ok){ + Log("FindTrackLengthInWater Tool: No MaxTotalHitsToDNN specified: assuming 1100, but this MUST match the value used for DNN training!",v_warning,verbosity); + maxhits0=1100; + } + Log("FindTrackLengthInWater Tool: max number of hits per event: "+to_string(maxhits0),v_debug,verbosity); + + // Get geometry variables from ANNIEEvent + // ============================= + get_ok = m_data->Stores.at("ANNIEEvent")->Header->Get("AnnieGeometry",anniegeom); + if(not get_ok){ + Log("FindTrackLengthInWater Tool: No Geometry in ANNIEEvent!",v_error,verbosity); + return false; + } + tank_radius = anniegeom->GetTankRadius()*100.; + std::cout<<"Tank radius is "<GetTankHalfheight()*100.; + std::cout<<"Tank halfheight is "<Branch("ievt", &ievt, "ievt/I"); - nu_eneNEW->Branch("neutrinoE", &trueNeuE, "neutrinoE/F"); - nu_eneNEW->Branch("trueKE", &trueE, "trueKE/F"); - nu_eneNEW->Branch("diffDirAbs2", &diffDirAbs2, "diffDirAbs2/F"); - nu_eneNEW->Branch("TrueTrackLengthInWater", &TrueTrackLengthInWater2, "TrueTrackLengthInWater/F"); - nu_eneNEW->Branch("TrueTrackLengthInMrd", &TrueTrackLengthInMrd2, "TrueTrackLengthInMrd/F"); - nu_eneNEW->Branch("recoDWallR2", &recoDWallR2, "recoDWallR2/F"); - nu_eneNEW->Branch("recoDWallZ2", &recoDWallZ2, "recoDWallZ2/F"); - nu_eneNEW->Branch("totalPMTs", &totalPMTs2, "totalPMTs/F"); - nu_eneNEW->Branch("totalLAPPDs", &totalLAPPDs2, "totalLAPPDs/F"); - nu_eneNEW->Branch("lambda_max_2", &lambda_max_2, "lambda_max_2/F"); - nu_eneNEW->Branch("dirX",&dirX2, "dirX/F"); - nu_eneNEW->Branch("dirY",&dirY2, "dirY/F"); - nu_eneNEW->Branch("dirZ",&dirZ2, "dirZ/F"); - nu_eneNEW->Branch("vtxX",&vtxX2, "vtxX/F"); - nu_eneNEW->Branch("vtxY",&vtxY2, "vtxY/F"); - nu_eneNEW->Branch("vtxZ",&vtxZ2, "vtxZ/F"); - - //----------- read the tree from file: - //deny_access=1; - Int_t run, event, nhits, trigger,recoStatus; - double vtxX,vtxY,vtxZ,dirX,dirY,dirZ,TrueTrackLengthInMrd,TrueTrackLengthInWater,TrueNeutrinoEnergy,trueEnergy,TrueMomentumTransfer,TrueMuonAngle; - std::string *TrueInteractionType = 0; - std::vector *digitX=0; std::vector *digitY=0; std::vector *digitZ=0; - std::vector *digitT=0; std::vector *digitType=0; - - regTree->GetEntry(currententry); - - regTree->SetBranchAddress("run", &run); - regTree->SetBranchAddress("event", &event); - regTree->SetBranchAddress("trueEnergy", &trueEnergy); - regTree->SetBranchAddress("TrueNeutrinoEnergy", &TrueNeutrinoEnergy); - regTree->SetBranchAddress("trigger", &trigger); - regTree->SetBranchAddress("nhits", &nhits); - regTree->SetBranchAddress("vtxX", &vtxX); - regTree->SetBranchAddress("vtxY", &vtxY); - regTree->SetBranchAddress("vtxZ", &vtxZ); - regTree->SetBranchAddress("dirX", &dirX); - regTree->SetBranchAddress("dirY", &dirY); - regTree->SetBranchAddress("dirZ", &dirZ); - regTree->SetBranchAddress("digitT", &digitT); - regTree->SetBranchAddress("digitX", &digitX); - regTree->SetBranchAddress("digitY", &digitY); - regTree->SetBranchAddress("digitZ", &digitZ); - regTree->SetBranchAddress("digitType", &digitType); - regTree->SetBranchAddress("recoStatus", &recoStatus); - regTree->SetBranchAddress("TrueInteractionType", &TrueInteractionType); - regTree->SetBranchAddress("TrueTrackLengthInMrd", &TrueTrackLengthInMrd); - regTree->SetBranchAddress("TrueTrackLengthInWater", &TrueTrackLengthInWater); - regTree->SetBranchAddress("TrueMomentumTransfer", &TrueMomentumTransfer); - regTree->SetBranchAddress("TrueMuonAngle", &TrueMuonAngle); - - double lambda_min = 10000000; double lambda_max = -99999999.9; double lambda = 0; - int totalPMTs=0; int totalLAPPDs=0; recoDWallR2=0; recoDWallZ2=0; diffDirAbs2=0; - double lambda_vec[1100]={0.}; double digitt[1100]={0.}; - - std::cout<<"currententry: "<0.){ - //std::cout<<"current entry: "<Stores.at("EnergyReco")->Delete();//clear the last entry from RAM + + // See if this event passes selection: we have several potential cuts. + // First check the EventCutstatus: event was in fiducial volume, had a stopping muon etc. + bool EventCutstatus; + get_ok = m_data->Stores.at("RecoEvent")->Get("EventCutStatus",EventCutstatus); + if(not get_ok){ + Log("FindTrackLengthInWater Tool: No EventCutStatus in the ANNIEEvent!",v_error,verbosity); + return false; + } + if(not EventCutstatus){ + Log("FindTrackLengthInWater Tool: Event did not pass the reconstruction selection cuts, skipping",v_message, verbosity); + return true; + } + count2++; + + // Check for the reconstructed tank vertex + RecoVertex* theExtendedVertex=nullptr; + get_ok = m_data->Stores.at("RecoEvent")->Get("ExtendedVertex", theExtendedVertex); + if((get_ok==0)||(theExtendedVertex==nullptr)){ + Log("FindTrackLengthInWater Tool: Failed to retrieve the ExtendedVertex from RecoEvent Store!",v_error,verbosity); + return false; + } + Int_t recoStatus = theExtendedVertex->GetStatus(); + double recoVtxFOM = theExtendedVertex->GetFOM(); + Log("FindTrackLengthInWater Tool: recoVtxFOM="+to_string(recoVtxFOM),v_debug,verbosity); + if(recoVtxFOM<=0.){ + Log("FindTrackLengthInWater Tool: Vertex reconstruction failed, skipping",v_message,verbosity); + return true; + } + Log("FindTrackLengthInWater Tool: Vertex reconstruction cuts passed",v_debug,verbosity); + count3++; + + //Get the reconstructed MRDTrackLength from CStore + std::vector> MrdTimeClusters; + std::vector* theMrdTracks; + int numtracksinev; + double MRDTrackLength=-9999; + Position StartVertex; + Position StopVertex; + //Load the mrd time clusters from CStore + bool get_clusters = m_data->CStore.Get("MrdTimeClusters",MrdTimeClusters); + if(!get_clusters){ + std::cout << "FindTrackLengthInWater tool: No MRD time clusters found. Will be no tracks." << std::endl; + return false; + } + + //loop through the time clusters and get the mrd tracks + for(int i=0; i < (int) MrdTimeClusters.size(); i++){ + m_data->Stores["MRDTracks"]->Get("MRDTracks",theMrdTracks); + m_data->Stores["MRDTracks"]->Get("NumMrdTracks",numtracksinev); + //loop through the mrd tracks and get the necessary info to calculate MRDTrackLength for each track + for(int tracki=0; trackiat(tracki)); + int TrackEventID = -1; + //get track properties that are needed for the through-going muon selection + thisTrackAsBoostStore->Get("MrdSubEventID",TrackEventID); + if(TrackEventID!= i) continue; + thisTrackAsBoostStore->Get("StartVertex",StartVertex); + thisTrackAsBoostStore->Get("StopVertex",StopVertex); + MRDTrackLength = sqrt(pow((StopVertex.X()-StartVertex.X()),2)+pow(StopVertex.Y()-StartVertex.Y(),2)+pow(StopVertex.Z()-StartVertex.Z(),2)) * 100.0; + } + } + //Don't use events without reconstructed track length in the MRD + if(MRDTrackLength<=0.){ + Log("FindTrackLengthInWater Tool: MRD reconstruction failed, skipping",v_message,verbosity); + return true; + } + count1++; + + // That's all the cuts! + // ==================== + //get info for mc vertex + RecoVertex *trueVtx = 0; + auto get_truevtx = m_data->Stores.at("RecoEvent")->Get("TrueVertex", trueVtx); + + //Get mc track length in the mrd stored in the RecoEvent store from previous tools + double TrueTrackLengthInMrd; + auto get_MRDTrackLength = m_data->Stores.at("RecoEvent")->Get("TrueTrackLengthInMRD", TrueTrackLengthInMrd); + Log("FindTrackLengthInWater Tool: TrueTrackLengthInMrd="+to_string(TrueTrackLengthInMrd),v_debug,verbosity); + + // proceed with getting event info we need to calculate the lambda values for the track length in water reconstruction + double vtxX,vtxY,vtxZ,dirX,dirY,dirZ,TrueNeutrinoEnergy,trueEnergy; + std::vector digitX; std::vector digitY; std::vector digitZ; + std::vector digitT; + + // get reconstructed vertex and direction info + Log("FindTrackLengthInWater Tool: Getting reco vertex and direction components",v_debug,verbosity); + vtxX = theExtendedVertex->GetPosition().X(); + vtxY = theExtendedVertex->GetPosition().Y(); + vtxZ = theExtendedVertex->GetPosition().Z(); + dirX = theExtendedVertex->GetDirection().X(); + dirY = theExtendedVertex->GetDirection().Y(); + dirZ = theExtendedVertex->GetDirection().Z(); + + // get additional primary muon info from RecoEvent store + Log("FindTrackLengthInWater Tool: Getting primary muon info",v_debug,verbosity); + double TrueTrackLengthInWater; + auto get_tankTrackLength = m_data->Stores.at("RecoEvent")->Get("TrueTrackLengthInWater", TrueTrackLengthInWater); + if(not get_tankTrackLength){ + Log("FindTrackLengthInWater Tool: Failed to retrieve the TrueTrackLengthInWater from RecoEvent Store!",v_error,verbosity); + return false; + } + auto get_muonMCEnergy = m_data->Stores.at("RecoEvent")->Get("TrueMuonEnergy", trueEnergy); + if(not get_muonMCEnergy){ + Log("FindTrackLengthInWater Tool: Failed to retrieve the TrueMuonEnergy from RecoEvent Store!",v_error,verbosity); + return false; + } + + /*// Get neutrino info needed for neutrino energy BDT training + get_ok = m_data->Stores.at("GenieEvent")->Get("TrueNeutrinoEnergy", TrueNeutrinoEnergy); + if(not get_ok){ + Log("FindTrackLengthInWater Tool: Failed to retrieve TrueNeutrinoEnergy!",v_error,verbosity); + return false; + } + */ + + // get digits from RecoDigit store + std::vector* digitList; + auto get_digits=m_data->Stores.at("RecoEvent")->Get("RecoDigit", digitList); + if(not get_digits){ + Log("FindTrackLengthInWater Tool: Failed to retrieve the RecoDigit from RecoEvent Store!",v_error,verbosity); + return false; + } + // Extract the PMT & LAPPD digit information + // =============================== + int totalPMTs =0; // number of digits from PMT hits in the event + int totalLAPPDs = 0; // number of digits from LAPPD hits in the event + //loop through all digits + for(RecoDigit &adigit : *digitList){ + digitX.push_back(adigit.GetPosition().X()); + digitY.push_back(adigit.GetPosition().Y()); + digitZ.push_back(adigit.GetPosition().Z()); + digitT.push_back(adigit.GetCalTime()); + if(adigit.GetDigitType()==0){totalPMTs+=1;} //when the digit type is zero we have a PMT digit + else{totalLAPPDs+=1;}// when it is 1 we have LAPPD + } + Log("FindTrackLengthInWater Tool: Got "+to_string(totalPMTs)+" PMT digits; "+to_string(digitT.size()) + +" total digits so far",v_debug,verbosity); + Log("FindTrackLengthInWater Tool: Got "+to_string(totalLAPPDs)+" LAPPD digits; "+to_string(digitT.size()) + +" total digits",v_debug,verbosity); + + // Estimate the track length in the tank + // ===================================== + //calculate diff dir with (0,0,1) double diffDirAbs0 = TMath::ACos(dirZ)*TMath::RadToDeg(); - //cout<<"diffDirAbs0: "<at(k)<<" | "<at(k)<at(k); - if( (digitType->at(k)) == "PMT8inch"){ totalPMTs++; } - if( (digitType->at(k)) == "lappd_v0"){ totalLAPPDs++; } - - //------ Find rack Length as the distance between the reconstructed vertex last photon emission point ----/ - lambda = find_lambda(vtxX,vtxY,vtxZ,dirX,dirY,dirZ,digitX->at(k),digitY->at(k),digitZ->at(k),42.); + float diffDirAbs2=diffDirAbs0/90.; + double recoVtxR2 = vtxX*vtxX + vtxZ*vtxZ; + double recoDWallR = tank_radius-TMath::Sqrt(recoVtxR2); + double recoDWallZ = tank_halfheight/2.-TMath::Abs(vtxY); + //Get true vtx and dir (optional) + double truevtxX=trueVtx->GetPosition().X(); + double truevtxY=trueVtx->GetPosition().Y(); + double truevtxZ=trueVtx->GetPosition().Z(); + double truedirX=trueVtx->GetDirection().X(); + double truedirY=trueVtx->GetDirection().Y(); + double truedirZ=trueVtx->GetDirection().Z(); + // Estimate the track length + // ========================= + Log("FindTrackLengthInWater Tool: Estimating track length in tank",v_debug,verbosity); + double lambda_min = 10000000; double lambda_max = -99999999.9; double lambda = 0; + std::vector lambda_vector; + for(int k=0; k= lambda_max ){ lambda_max = lambda; } - lambda_vec[k]=lambda; - //m_data->Stores["ANNIEEvent"]->Set("WaterRecoTrackLength",lambda_max); + lambda_vector.push_back(lambda); } - //std::cout<<"the track length in the water tank (1st approx) is: "<Fill(); - } - } - currententry++; - - if(currententry==regTree->GetEntries()) m_data->vars.Set("StopLoop",1); - //if(currententry==5) m_data->vars.Set("StopLoop",1); + + // Post-processing of variables to store + // ===================================== + float recoDWallR2 = recoDWallR/tank_radius; + float recoDWallZ2 = recoDWallZ/tank_halfheight*2.; + float TrueTrackLengthInWater2 = TrueTrackLengthInWater*100.; // convert to [cm] + float TrueTrackLengthInMrd2 = TrueTrackLengthInMrd; // it is already in [cm] + // we need to normalise the digit time and lambda vectors to fixed dimensions to match the MaxTotalHitsToDNN + + lambda_vector.resize(maxhits0); + + digitT.resize(maxhits0); + + //put MaxTotalHitsToDNN in store for use in the next tools + m_data->Stores.at("EnergyReco")->Set("MaxTotalHitsToDNN",maxhits0); + + // put the last successfully processed event number in the EnergyReco store as well. + // write the current event to file + uint32_t EventNumber; + get_ok = m_data->Stores.at("ANNIEEvent")->Get("EventNumber", EventNumber); + + + // Put these variables in the EnergyReco BoostStore + // ================================================ + Log("FindTrackLengthInWater Tool: putting event "+to_string(EventNumber)+" into the EnergyReco store",v_debug,verbosity); + m_data->Stores.at("EnergyReco")->Set("ThisEvtNum",EventNumber); + //std::cout<<"This is the Eventnumber you are looking for"<Stores.at("EnergyReco")->Set("lambda_vec",lambda_vector); + m_data->Stores.at("EnergyReco")->Set("digit_ts_vec",digitT); + m_data->Stores.at("EnergyReco")->Set("lambda_max",lambda_max); + m_data->Stores.at("EnergyReco")->Set("num_pmt_hits",totalPMTs); + m_data->Stores.at("EnergyReco")->Set("num_lappd_hits",totalLAPPDs); + m_data->Stores.at("EnergyReco")->Set("TrueTrackLengthInWater",TrueTrackLengthInWater2); + //m_data->Stores.at("EnergyReco")->Set("trueNeuE",TrueNeutrinoEnergy); + m_data->Stores.at("EnergyReco")->Set("trueMuonEnergy",trueEnergy); + m_data->Stores.at("EnergyReco")->Set("diffDirAbs",diffDirAbs2); + //m_data->Stores.at("EnergyReco")->Set("TrueTrackLengthInMrd",TrueTrackLengthInMrd2); + m_data->Stores.at("EnergyReco")->Set("recoDWallR",recoDWallR2); + m_data->Stores.at("EnergyReco")->Set("recoDWallZ",recoDWallZ2); + m_data->Stores.at("EnergyReco")->Set("dirVec",theExtendedVertex->GetDirection()); + m_data->Stores.at("EnergyReco")->Set("vtxVec",theExtendedVertex->GetPosition()); + m_data->Stores.at("EnergyReco")->Set("recoTrackLengthInMrd",MRDTrackLength); + + if(not fDoTraining){ + Log("FindTrackLengthInWater Tool: Not doing training, so BoostStore will not be saved locally",v_error,verbosity); + } +// only save all data to disk when training +if(fDoTraining){ + m_data->Stores.at("EnergyReco")->Save("EnergyReco.bs"); // write this Energy Reco entry to disk. If the file exists, it appends as a new entry. + } +//---------------------------------------------------------------------------------------------------------------------------------------------------------// + if(lambda_vector.size()!=maxhits0){ + Log("FindTrackLengthInWater Tool: Error! lambdavector size is not maxhits0! Check dimensions!",v_error,verbosity); + return false; + } + if(digitT.size()!=maxhits0){ + Log("FindTrackLengthInWater Tool: Error! digitT size is not maxhits0! Check dimensions!",v_error,verbosity); + return false; + } + + // Write to .csv file + // ================== + if(not csvfile.is_open()){ + Log("FindTrackLengthInWater Tool: output file is closed, skipping write",v_debug,verbosity); + return true; + } + + for(int i=0; iWrite(); + m_data->Stores.at("EnergyReco")->Close(); + if(csvfile.is_open()) csvfile.close(); + Log("FindTrackLengthInWater Tool: processed "+to_string(count1)+" events",v_message,verbosity); + std::cout<<"processed a total of "< #include -#include #include "ANNIEalgorithms.h" #include "Tool.h" @@ -25,25 +24,24 @@ class FindTrackLengthInWater: public Tool { private: - std::string infile; - TFile* file; - TTree* regTree; - TTree * nu_eneNEW; - ExampleRoot* Data; - int maxhits0=1100; - long currententry; - long NumEvents; bool first=1; bool deny_access=0; - double diffDirAbs2=0; double diffDirAbs=0; - double recoDWallR2=0; double recoDWallZ2=0; - int count1=0; - + // counters to keep track of cut efficiencies + int count1=0, count2=0, count3=0, count4=0; std::ofstream csvfile; - std::string myfile; - std::string outputdir=""; - bool writefile=false; - TFile* outputFile; + Geometry* anniegeom=nullptr; + double tank_radius; + double tank_halfheight; + int fDoTraining=0; + + // verbosity levels: if 'verbosity' < this level, the message type will be logged. + int verbosity=1; + int v_error=0; + int v_warning=1; + int v_message=2; + int v_debug=3; + std::string logmessage; + int get_ok; }; diff --git a/UserTools/FindTrackLengthInWater/README.md b/UserTools/FindTrackLengthInWater/README.md index 915425343..4fdd4ff5e 100644 --- a/UserTools/FindTrackLengthInWater/README.md +++ b/UserTools/FindTrackLengthInWater/README.md @@ -1,7 +1,72 @@ # FindTrackLengthInWater -FindTrackLengthInWater +The `FindTrackLengthInWater` tool uses the `find_lambda()` function to calculate the length between the reconstructed Vertex and the photon production Point(lambda) for each digit. These lengths are going to be used to reconstruct the track in the water in the next tool. In order to do that the tool requires information about the `ExtendedVertex` and the `RecoDigit` which are stored in the `RecoEvent` store from previous tools. Then the required information for the track length and energy reconstruction are passed to the next tools through the `EnergyReco` boost store. Also cuts are applied to exclude events with bad vertex reconstruction(recovtxFOM<=0) and with no reconstructed tracks in the MRD. -* Takes the .root file from reconstruction (e.g. reco26_5LAPPDs+128PMTs_extv3_9.root), calculates the track length as the distance between the first and last Cherenkov photon emission point along the track and stores variables for track length reconstruction using DNN in a .csv file. -* Set input/output file names at: configfiles/FindTrackLengthInWater/LoadRecoInputFile +## Data +### Input + +The following variables are obtained from the `ANNIEEvent` store: + +**EventNumber** `uint32_t` + +**AnnieGeometry**(Header) `Geometry` Retrieve the tank radius and the tank halfheight from the geometry object of the `ANNIEEvent` store + +The following variables are obtained from the `RecoEvent` store: + +**ExtendedVertex** `RecoVertex` Contains the reconstructed Vertex information + +**MrdTimeClusters** `vector>` One vector for each subevent containing the digit IDs for clustered hits in each subevent + +**TrueTrackLengthInWater** `double` MC track length in the water is retrieved from the `RecoEvent` store + +**TrueMuonEnergy** `double` MC muon energy is retrieved from the `RecoEvent` store + +**RecoDigit** `std::vector` One vector for each event containing the RecoDigit objects for each event + +### Output + +The following variables are passed on to the next tool via the `EnergyReco` store: + +**MaxTotalHitsToDNN** `int` + +**ThisEvtNum** `uint32_t` + +**lambda_vec** `std::vector` Vector with all the lambda values for each event + +**digit_ts_vec** `std::vector` Vector with the time of all the digits of each event + +**lambda_max** `double` The distance between the reconstructed vertex and last Cherenkov photon emission point along the track + +**num_pmt_hits** `int` Total number of pmt digits + +**num_lappd_hits** `int` Total number of lappd digits + +**TrueTrackLengthInWater** `float` MC track length in the water + +**trueMuonEnergy** `double` MC muon energy + +**diffDirAbs** `double` Angle difference between the reconstructed z direction and the beam direction at (0,0,1) + +**recoDWallR** `double` Radial distance of the reconstructed vertex from the walls of the tank + +**recoDWallZ** `double` Axial distance of the reconstructed vertex from the walls of the tank + +**dirVec** `Direction` Direction of the reconstructed vertex + +**vtxVec** `Position` Position of the reconstructed vertex + +**recoTrackLengthInMrd** Track length of a reconstructed track in the MRD found by the `FindMrdTracks` tool + +## Configuration + +``` +MaxTotalHitsToDNN 1100 +This is set to 1100 as a reasonable upper limit to the total number of hits that might happen in a single event in order to have a fixed number of features for the DNN + +OutputDataFile data_for_trackLength_training.csv +There exists the capabillity to create an output .csv file with the necessary information for the DNN for any inpdependent analysis or testing + +DoTraining bool +It is by default set to 0. Only when you want to run the training toolchain you need to set it to 1 +``` diff --git a/UserTools/PlotsTrackLengthAndEnergy/PlotsTrackLengthAndEnergy.cpp b/UserTools/PlotsTrackLengthAndEnergy/PlotsTrackLengthAndEnergy.cpp new file mode 100644 index 000000000..a8b33ab95 --- /dev/null +++ b/UserTools/PlotsTrackLengthAndEnergy/PlotsTrackLengthAndEnergy.cpp @@ -0,0 +1,121 @@ +#include "PlotsTrackLengthAndEnergy.h" +#include "TCanvas.h" +#include "TMath.h" +#include "TH2F.h" +#include "TLegend.h" +#include "TAxis.h" +#include "TLine.h" + +PlotsTrackLengthAndEnergy::PlotsTrackLengthAndEnergy():Tool(){} + + +bool PlotsTrackLengthAndEnergy::Initialise(std::string configfile, DataModel &data){ + + /////////////////// Useful header /////////////////////// + if(configfile!="") m_variables.Initialise(configfile); // loading config file + //m_variables.Print(); + + m_data= &data; //assigning transient data pointer + ///////////////////////////////////////////////////////////////// + + return true; +} + + +bool PlotsTrackLengthAndEnergy::Execute(){ + + BoostStore EnergyReco(true,2); + + EnergyReco.Initialise("EnergyRecoStore.bs"); + + unsigned long n_entries; + bool get_ok = EnergyReco.Header->Get("TotalEntries", n_entries); + if(not get_ok) { + Log("PlotsTrackLengthAndEnergy Tool: EnergyRecoStore file does not exist, run the EnergyRecoPredict toolchain first!",v_error,verbosity); + return false; + } + std::cout<<"got total entries; "< +#include + +#include "Tool.h" + + +/** + * \class PlotsTrackLengthAndEnergy + * + * This is a blank template for a Tool used by the script to generate a new custom tool. Please fill out the description and author information. +* +* $Author: B.Richards $ +* $Date: 2019/05/28 10:44:00 $ +* Contact: b.richards@qmul.ac.uk +*/ +class PlotsTrackLengthAndEnergy: public Tool { + + + public: + + PlotsTrackLengthAndEnergy(); ///< Simple constructor + bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools. + bool Execute(); ///< Execute function used to perform Tool purpose. + bool Finalise(); ///< Finalise function used to clean up resources. + + + private: + +// verbosity levels: if 'verbosity' < this level, the message type will be logged. + int verbosity=1; + int v_error=0; + int v_warning=1; + int v_message=2; + int v_debug=3; + std::string logmessage; + + + +}; + + +#endif diff --git a/UserTools/PlotsTrackLengthAndEnergy/README.md b/UserTools/PlotsTrackLengthAndEnergy/README.md new file mode 100644 index 000000000..c3f1db6f5 --- /dev/null +++ b/UserTools/PlotsTrackLengthAndEnergy/README.md @@ -0,0 +1,23 @@ +# PlotsTrackLengthAndEnergy + +The `PlotsTrackLengthAndEnergy` tool is used to plot the reconstructed track length and energy from the DNN and the BDT. Therefore this tool should be run only after running the `EnergyRecoPredict` toolchain. This tool loads the `EnergyReco` boost store that is being saved locally by the `BDTMuonEnergyPredict` tool in order to make the plots we want. + +## Data + +### Input + +The following variables are obtained from the `EnergyRecoStore`: + +**lambda_max** `double` The distance between the reconstructed vertex and last Cherenkov photon emission point along the track + +**TrueTrackLengthInWater** `float` MC track length in the water + +**DNNRecoLength** `double` The reconstructed track length in water for a single event + +**trueMuonEnergy** `double` MC muon energy + +**BDTMuonEnergy** `double` The reconstructed energy of the muon for a single event + +## Configuration + +There aren't any configuration variables for PlotsTrackLengthAndEnergy. diff --git a/UserTools/Unity.h b/UserTools/Unity.h index db441a61d..dc6e501b7 100644 --- a/UserTools/Unity.h +++ b/UserTools/Unity.h @@ -167,5 +167,6 @@ #include "FilterEvents.h" #include "Stage1DataBuilder.h" #include "BeamFetcherV2.h" +#include "PlotsTrackLengthAndEnergy.h" #include "SaveConfigInfo.h" #include "ReadConfigInfo.h" diff --git a/configfiles/EnergyReco/DummyToolConfig b/configfiles/EnergyReco/DummyToolConfig deleted file mode 100644 index 95cad88d2..000000000 --- a/configfiles/EnergyReco/DummyToolConfig +++ /dev/null @@ -1,3 +0,0 @@ -# Dummy config file - -verbose 2 \ No newline at end of file diff --git a/configfiles/EnergyReco/Erecoconfig b/configfiles/EnergyReco/Erecoconfig deleted file mode 100644 index ae289ebdf..000000000 --- a/configfiles/EnergyReco/Erecoconfig +++ /dev/null @@ -1,6 +0,0 @@ -#PythonScript BDT_MuonEnergyReco -PythonScript BDT_NeutrinoEnergyReco - -InitialiseFunction Initialise -ExecuteFunction Execute -FinaliseFunction Finalise diff --git a/configfiles/EnergyReco/Predict/BDTMuonEnergyPredictConfig b/configfiles/EnergyReco/Predict/BDTMuonEnergyPredictConfig new file mode 100644 index 000000000..6eb0a32e5 --- /dev/null +++ b/configfiles/EnergyReco/Predict/BDTMuonEnergyPredictConfig @@ -0,0 +1,9 @@ +# This config file stores variables that may be retrieved by the python script +# Access them by retrieving the corresponding variable from the 'Config' store + +PythonScript BDTMuonEnergyPredict + +verbose 1 + +# BDT model file to load +BDTMuonModelFile finalized_BDTmodel_forMuonEnergy.sav diff --git a/configfiles/EnergyReco/Predict/ClusterFinderConfig b/configfiles/EnergyReco/Predict/ClusterFinderConfig new file mode 100644 index 000000000..05600b7fb --- /dev/null +++ b/configfiles/EnergyReco/Predict/ClusterFinderConfig @@ -0,0 +1,12 @@ +# ClusterFinder Config File + +HitStore MCHits #Either MCHits or Hits (accessed in ANNIEEvent store) +OutputFile Run1415S0_AllPMTs_ClusterFinder #Output root prefix name for the current run +ClusterFindingWindow 50 # in ns, size of the window used to "clusterize" +AcqTimeWindow 4000 # in ns, size of the acquisition window +ClusterIntegrationWindow 50 # in ns, all hits with +/- 1/2 of this window are considered in the cluster +MinHitsPerCluster 10 # group of hits are considered clusters above this amount of hits + +verbosity 0 #verbosity of the application +end_of_window_time_cut 0.95 # from 0 to 1, length of the window you want to loop over with respect to acq. window (1 for full window, 0.95 for 95% from the start) +Plots2D 0 #2D charge-vs-time plot to be drawn? diff --git a/configfiles/EnergyReco/Predict/DNNTrackLengthPredictConfig b/configfiles/EnergyReco/Predict/DNNTrackLengthPredictConfig new file mode 100644 index 000000000..7dd7ac486 --- /dev/null +++ b/configfiles/EnergyReco/Predict/DNNTrackLengthPredictConfig @@ -0,0 +1,13 @@ +# This config file stores variables that may be retrieved by the python script +# Also this is a caller config file + +PythonScript DNNTrackLengthPredict + +verbose 1 + +# DNN model, weights etc +TrackLengthWeightsFile UserTools/DNNTrackLength/stand_alone/weights/weights_bets.hdf5 + +# Scaling Parameters + +ScalingVarsBoostStoreFile Data_Energy_Reco/ScalingVarsStore.bs diff --git a/configfiles/EnergyReco/Predict/DigitBuilderConfig b/configfiles/EnergyReco/Predict/DigitBuilderConfig new file mode 100644 index 000000000..8eed149e9 --- /dev/null +++ b/configfiles/EnergyReco/Predict/DigitBuilderConfig @@ -0,0 +1,11 @@ +# DigitBuilder config file + +verbosity 0 +ParametricModel 1 +#Reading in MC files +IsMC 1 +DigitChargeThr 20 +# There are three configurations: "PMT_only", "LAPPD_only", "All" +PhotoDetectorConfiguration All +#File must be in /pnfs/ space when loading on Fermilab cluster +LAPPDIDFile configfiles/EnergyReco/Predict/LAPPDIDs.txt diff --git a/configfiles/EnergyReco/Predict/EventSelectorConfig b/configfiles/EnergyReco/Predict/EventSelectorConfig new file mode 100644 index 000000000..ca6032811 --- /dev/null +++ b/configfiles/EnergyReco/Predict/EventSelectorConfig @@ -0,0 +1,15 @@ +# EventSelector config file + +verbosity 0 +MCPMTVolCut 0 +MCFVCut 1 +MCMRDCut 1 +MCPiKCut 0 +MCIsMuonCut 1 +MRDRecoCut 0 +RecoPMTVolCut 0 +RecoFVCut 0 +NHitCut 1 +PromptTrigOnly 1 +SaveStatusToStore 1 +Triggerword -1 diff --git a/configfiles/EnergyReco/Predict/FindMrdTracksConfig b/configfiles/EnergyReco/Predict/FindMrdTracksConfig new file mode 100644 index 000000000..7a9d11f7d --- /dev/null +++ b/configfiles/EnergyReco/Predict/FindMrdTracksConfig @@ -0,0 +1,12 @@ +# FindMrdTracks Config File +# all variables retrieved with m_variables.Get() must be defined here! + +verbosity 0 +IsData 0 +OutputDirectory . +OutputFile STEC_MRDTracks_cluster40ns +DrawTruthTracks 0 # whether to add MC Truth track info for drawing in MrdPaddlePlot Tool + ## note you need to run that tool to actually view the tracks! +WriteTracksToFile 0 # should the track information be written to a ROOT-file? +SelectTriggerType 0 #should the loaded data be filtered by trigger type? +TriggerType Beam #options: Cosmic, Beam, No Loopback diff --git a/configfiles/EnergyReco/Predict/FindTrackLengthInWaterConfig b/configfiles/EnergyReco/Predict/FindTrackLengthInWaterConfig new file mode 100644 index 000000000..2322e20c5 --- /dev/null +++ b/configfiles/EnergyReco/Predict/FindTrackLengthInWaterConfig @@ -0,0 +1,5 @@ +#FindTrackLengthInWater Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbosity 1 +MaxTotalHitsToDNN 1100 +OutputDataFile beamlikeEvts.csv diff --git a/configfiles/EnergyReco/Predict/HitCleanerConfig b/configfiles/EnergyReco/Predict/HitCleanerConfig new file mode 100644 index 000000000..36d7f0137 --- /dev/null +++ b/configfiles/EnergyReco/Predict/HitCleanerConfig @@ -0,0 +1,20 @@ +# HitCleaner config file + +verbosity 5 +Config 2 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +PmtMinPulseHeight 20 #minimum pulse height +PmtNeighbourRadius 60 #digit neighbouring distance [cm] +PmtMinNeighbourDigits 2 #minimum neighbour digits +PmtClusterRadius 60 #digit clustering distance [cm] +PmtTimeWindowN 10 #neighbouring time window [ns] +PmtTimeWindowC 10 #clustering time window [ns] +PmtMinHitsPerCluster 4 #number of hits per cluster +LappdMinPulseHeight 0 #minimum pulse height +LappdNeighbourRadius 25 #digit neighbouring distance [cm] +LappdMinNeighbourDigits 20 #minimum neighbour digits +LappdClusterRadius 25 #digit clustering distance [cm] +LappdTimeWindowN 1 #neighbouring time window [ns] +LappdTimeWindowC 1 #clustering time window [ns] +LappdMinHitsPerCluster 20 #number of digits per cluster +#MinClusterDigits 50 #minimum clustered digits (LAPPD+PMT) +MinClusterDigits 10 #minimum clustered digits (PMT-only) diff --git a/configfiles/EnergyReco/Predict/LAPPDIDs.txt b/configfiles/EnergyReco/Predict/LAPPDIDs.txt new file mode 100644 index 000000000..776b28912 --- /dev/null +++ b/configfiles/EnergyReco/Predict/LAPPDIDs.txt @@ -0,0 +1,5 @@ +11 +13 +14 +15 +17 diff --git a/configfiles/EnergyReco/Predict/LoadWCSimConfig b/configfiles/EnergyReco/Predict/LoadWCSimConfig new file mode 100644 index 000000000..001d2e929 --- /dev/null +++ b/configfiles/EnergyReco/Predict/LoadWCSimConfig @@ -0,0 +1,15 @@ +#LoadWCSim Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbose 0 + +InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standard/tank/pmt/wcsim_0* #James's files + +TriggerType Beam +SplitSubTriggers 1 + +WCSimVersion 10 ## should reflect the WCSim version of the files being loaded +HistoricTriggeroffset 0 ## time offset of digits relative to the trigger +#UseDigitSmearedTime 1 ## whether to use smeared digit time (T), or true time of first photon (F) +LappdNumStrips 60 ## num channels to construct from each LAPPD +LappdStripLength 100 ## relative x position of each LAPPD strip, for dual-sided readout [mm] +LappdStripSeparation 10 ## stripline separation, for calculating relative y position of each LAPPD strip [mm] diff --git a/configfiles/EnergyReco/Predict/LoadWCSimLAPPDConfig b/configfiles/EnergyReco/Predict/LoadWCSimLAPPDConfig new file mode 100644 index 000000000..2b82ff922 --- /dev/null +++ b/configfiles/EnergyReco/Predict/LoadWCSimLAPPDConfig @@ -0,0 +1,10 @@ +#LoadWCSimLAPPD Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbose 0 + +InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standard/tank/lappd/wcsim_lappd_0* #James's files + +WCSimVersion 10 ## should reflect the WCSim version of the files being loaded +InnerStructureRadius 1.3545 ## octagonal inner structure radius in m (from drawings 106.64") +HistoricTriggeroffset 0 +DrawDebugGraphs 0 ## whether to draw TPolyMarker3D's of hits diff --git a/configfiles/EnergyReco/Predict/MCParticlePropertiesConfig b/configfiles/EnergyReco/Predict/MCParticlePropertiesConfig new file mode 100644 index 000000000..752200268 --- /dev/null +++ b/configfiles/EnergyReco/Predict/MCParticlePropertiesConfig @@ -0,0 +1 @@ +verbosity 0 diff --git a/configfiles/EnergyReco/Predict/MCRecoEventLoaderConfig b/configfiles/EnergyReco/Predict/MCRecoEventLoaderConfig new file mode 100644 index 000000000..193bfd196 --- /dev/null +++ b/configfiles/EnergyReco/Predict/MCRecoEventLoaderConfig @@ -0,0 +1,12 @@ +# MCRecoEventLoader config file + +verbosity 0 +#Get Pion/Kaon counts from MC Truth +GetPionKaonInfo 1 +#Particle ID you want to load from parent particles; 13 for muon +ParticleID 13 +#coordinate shifts needed to put particles in tank origin coordinates +xshift 0.0 +yshift 14.46469 +zshift -168.1 + diff --git a/configfiles/EnergyReco/Predict/PlotsTrackLengthAndEnergyConfig b/configfiles/EnergyReco/Predict/PlotsTrackLengthAndEnergyConfig new file mode 100644 index 000000000..ef7a3d415 --- /dev/null +++ b/configfiles/EnergyReco/Predict/PlotsTrackLengthAndEnergyConfig @@ -0,0 +1 @@ +#PlotsTrackLengthAndEnergy config file diff --git a/configfiles/EnergyReco/Predict/TimeClusteringConfig b/configfiles/EnergyReco/Predict/TimeClusteringConfig new file mode 100644 index 000000000..03419833f --- /dev/null +++ b/configfiles/EnergyReco/Predict/TimeClusteringConfig @@ -0,0 +1,12 @@ +#TimeClustering config file + +verbosity 1 +MinDigitsForTrack 4 +MaxMrdSubEventDuration 40 +MinSubeventTimeSep 40 +MakeMrdDigitTimePlot 0 +LaunchTApplication 0 +IsData 0 +#OutputROOTFile TimeClustering_MRDTest28_cluster40ns +OutputROOTFile STEC_TimeClusteringOut +MapChankey_WCSimID ./configfiles/SimpleTankEnergyCalibrator/MRD_Chankey_WCSimID.dat diff --git a/configfiles/EnergyReco/ToolChainConfig b/configfiles/EnergyReco/Predict/ToolChainConfig similarity index 71% rename from configfiles/EnergyReco/ToolChainConfig rename to configfiles/EnergyReco/Predict/ToolChainConfig index db437bc5c..ad7c94761 100644 --- a/configfiles/EnergyReco/ToolChainConfig +++ b/configfiles/EnergyReco/Predict/ToolChainConfig @@ -15,9 +15,9 @@ service_publish_sec -1 service_kick_sec -1 ##### Tools To Add ##### -Tools_File configfiles/EnergyReco/ToolsConfig ## list of tools to run and their config files +Tools_File configfiles/EnergyReco/Predict/ToolsConfig ## list of tools to run and their config files ##### Run Type ##### -Inline 1 #-1 ## number of Execute steps in program, -1 infinite loop that is ended by user -Interactive 1 #0 ## set to 1 if you want to run the code interactively +Inline -1 ## number of Execute steps in program, -1 infinite loop that is ended by user +Interactive 0 ## set to 1 if you want to run the code interactively diff --git a/configfiles/EnergyReco/Predict/ToolsConfig b/configfiles/EnergyReco/Predict/ToolsConfig new file mode 100644 index 000000000..5926759d1 --- /dev/null +++ b/configfiles/EnergyReco/Predict/ToolsConfig @@ -0,0 +1,19 @@ +# predict machine learning algorithm +myLoadWCSim LoadWCSim configfiles/EnergyReco/Predict/LoadWCSimConfig +myLoadWCSimLAPPD LoadWCSimLAPPD configfiles/EnergyReco/Predict/LoadWCSimLAPPDConfig +myMCParticleProperties MCParticleProperties configfiles/EnergyReco/Predict/MCParticlePropertiesConfig +MCRecoEventLoader MCRecoEventLoader configfiles/EnergyReco/Predict/MCRecoEventLoaderConfig +DigitBuilder DigitBuilder configfiles/EnergyReco/Predict/DigitBuilderConfig +HitCleaner HitCleaner configfiles/EnergyReco/Predict/HitCleanerConfig +ClusterFinder ClusterFinder configfiles/EnergyReco/Predict/ClusterFinderConfig +TimeClustering TimeClustering configfiles/EnergyReco/Predict/TimeClusteringConfig +EventSelector EventSelector configfiles/EnergyReco/Predict/EventSelectorConfig +FindMrdTracks FindMrdTracks configfiles/EnergyReco/Predict/FindMrdTracksConfig +VtxSeedGenerator VtxSeedGenerator configfiles/EnergyReco/Predict/VtxSeedGeneratorConfig +VtxSeedFineGrid VtxSeedFineGrid configfiles/EnergyReco/Predict/VtxSeedFineGridConfig +VtxExtendedVertexFinder VtxExtendedVertexFinder configfiles/EnergyReco/Predict/VtxExtendedVertexFinderConfig + +FindTrackLengthInWater FindTrackLengthInWater configfiles/EnergyReco/Predict/FindTrackLengthInWaterConfig +DNNTrackLengthPredict PythonScript configfiles/EnergyReco/Predict/DNNTrackLengthPredictConfig +BDTMuonEnergyPredict PythonScript configfiles/EnergyReco/Predict/BDTMuonEnergyPredictConfig +#PlotsTrackLengthAndEnergy PlotsTrackLengthAndEnergy configfiles/EnergyReco/Predict/PlotsTrackLengthAndEnergyConfig diff --git a/configfiles/EnergyReco/Predict/VtxExtendedVertexFinderConfig b/configfiles/EnergyReco/Predict/VtxExtendedVertexFinderConfig new file mode 100644 index 000000000..59f74e5eb --- /dev/null +++ b/configfiles/EnergyReco/Predict/VtxExtendedVertexFinderConfig @@ -0,0 +1,9 @@ +# VtxPointPositionFinder config file + +verbosity 5 + +# Only one of the options below should be set to 1. All the others should be set to 0 +UseTrueVertexAsSeed 0 +UsePointVertexAsSeed 0 +FitAllOnSeedGrid 1 +UseMeanTimeAsSeed 0 diff --git a/configfiles/EnergyReco/Predict/VtxSeedFineGridConfig b/configfiles/EnergyReco/Predict/VtxSeedFineGridConfig new file mode 100644 index 000000000..674e643a0 --- /dev/null +++ b/configfiles/EnergyReco/Predict/VtxSeedFineGridConfig @@ -0,0 +1,7 @@ +verbosity 5 +useDirectionGrid 0 +useSimpleDirection 1 +useTrueDir 0 +useMRDTrack 0 +usePastResolution 0 +multiGrid 1 diff --git a/configfiles/EnergyReco/Predict/VtxSeedGeneratorConfig b/configfiles/EnergyReco/Predict/VtxSeedGeneratorConfig new file mode 100644 index 000000000..fe2688cc3 --- /dev/null +++ b/configfiles/EnergyReco/Predict/VtxSeedGeneratorConfig @@ -0,0 +1,7 @@ +# VtxSeedGenerator config file + +verbosity 3 +NumberOfSeeds 300 +UseSeedGrid 1 +SeedType 2 +UseFineGrid 0 diff --git a/configfiles/EnergyReco/README.md b/configfiles/EnergyReco/README.md index 5afa52c11..c7082bf29 100644 --- a/configfiles/EnergyReco/README.md +++ b/configfiles/EnergyReco/README.md @@ -1,25 +1,60 @@ # Configure files *********************** -#Description +# Description ********************** -Configure files are simple text files for passing variables to the Tools. - -Text files are read by the Store class (src/Store) and automatically asigned to an internal map for the relavent Tool to use. +The `EnergyReco` toolchains are used to reconstruct the track length in the water tank and the tank energy of a muon by using a Deep Learning Neural Network(DNN) and a Boosted Decision Tree(BDT) with Gradient Boost. The `Train_Test` toolchain is used to train and test the models that we want to use later for the reconstruction. The `Predict` toolchain is used after we've trained the models in order to reconstruct the track length in the water tank and the tank energy of a muon for different events. +Both of these toolchains require information from the digit, the MRD and the vertex reconstruction so the appropriate tools need to be used. ************************ -#Useage +# Tools ************************ -Any line starting with a "#" will be ignored by the Store, as will blank lines. - -Variables should be stored one per line as follows: - +************************ +# Train_Test +************************ -Name Value #Comments +The tools in the `Train_Test` toolchain are the following: + +* LoadWCSim +* LoadWCSimLAPPD +* MCParticleProperties +* MCRecoEventLoader +* DigitBuilder +* HitCleaner +* ClusterFinder +* TimeClustering +* EventSelector +* FindMrdTracks +* VtxSeedGenerator +* VtxSeedFineGrid +* VtxExtendedVertexFinder +* FindTrackLengthInWater +* DNNTrackLengthTrain_Test +* BDTMuonEnergyTrain_Test +************************ +# Predict +************************ -Note: Only one value is permitted per name and they are stored in a string stream and templated cast back to the type given. +The tools in the `Predict` toolchain are the following: + +* LoadWCSim +* LoadWCSimLAPPD +* MCParticleProperties +* MCRecoEventLoader +* DigitBuilder +* HitCleaner +* ClusterFinder +* TimeClustering +* EventSelector +* FindMrdTracks +* VtxSeedGenerator +* VtxSeedFineGrid +* VtxExtendedVertexFinder +* FindTrackLengthInWater +* DNNTrackLengthPredict +* BDTMuonEnergyPredict diff --git a/configfiles/EnergyReco/ToolsConfig b/configfiles/EnergyReco/ToolsConfig deleted file mode 100644 index 5bf9f59c1..000000000 --- a/configfiles/EnergyReco/ToolsConfig +++ /dev/null @@ -1 +0,0 @@ -EnergyReco PythonScript configfiles/EnergyReco/Erecoconfig diff --git a/configfiles/EnergyReco/Train_Test/BDTMuonEnergyTrainConfig b/configfiles/EnergyReco/Train_Test/BDTMuonEnergyTrainConfig new file mode 100644 index 000000000..d5fdbff87 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/BDTMuonEnergyTrainConfig @@ -0,0 +1,15 @@ +# This config file stores variables that may be retrieved by the python script +# Access them by retrieving the corresponding variable from the 'Config' store + +PythonScript BDTMuonEnergyTrain_Test + +verbosity 1 + +InitialiseFunction Initialise +ExecuteFunction Finalise +FinaliseFunction Execute + +# BDT model file to write +BDTMuonEnergyWeightsFile finalized_BDTmodel_forMuonEnergy.sav + +BDTMuonEnergyTrainingInputBoostStoreFile EnergyReco.bs diff --git a/configfiles/EnergyReco/Train_Test/ClusterFinderConfig b/configfiles/EnergyReco/Train_Test/ClusterFinderConfig new file mode 100644 index 000000000..05600b7fb --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/ClusterFinderConfig @@ -0,0 +1,12 @@ +# ClusterFinder Config File + +HitStore MCHits #Either MCHits or Hits (accessed in ANNIEEvent store) +OutputFile Run1415S0_AllPMTs_ClusterFinder #Output root prefix name for the current run +ClusterFindingWindow 50 # in ns, size of the window used to "clusterize" +AcqTimeWindow 4000 # in ns, size of the acquisition window +ClusterIntegrationWindow 50 # in ns, all hits with +/- 1/2 of this window are considered in the cluster +MinHitsPerCluster 10 # group of hits are considered clusters above this amount of hits + +verbosity 0 #verbosity of the application +end_of_window_time_cut 0.95 # from 0 to 1, length of the window you want to loop over with respect to acq. window (1 for full window, 0.95 for 95% from the start) +Plots2D 0 #2D charge-vs-time plot to be drawn? diff --git a/configfiles/EnergyReco/Train_Test/DNNTrackLengthTrainConfig b/configfiles/EnergyReco/Train_Test/DNNTrackLengthTrainConfig new file mode 100644 index 000000000..185227cd2 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/DNNTrackLengthTrainConfig @@ -0,0 +1,18 @@ +# This config file stores variables that may be retrieved by the python script +# Access them by retrieving the corresponding variable from the 'Config' store + +PythonScript DNNTrackLengthTrain_Test + +verbosity 1 + +InitialiseFunction Initialise +ExecuteFunction Finalise +FinaliseFunction Execute + +# DNN model, weights etc +TrackLengthOutputWeightsFile UserTools/DNNTrackLength/stand_alone/weights/weights_bets.hdf5 + +# Training data +TrackLengthTrainingInputBoostStoreFile EnergyReco.bs + +ScalingVarsBoostStoreFile Data_Energy_Reco/ScalingVarsStore.bs diff --git a/configfiles/EnergyReco/Train_Test/DigitBuilderConfig b/configfiles/EnergyReco/Train_Test/DigitBuilderConfig new file mode 100644 index 000000000..1dcc890b6 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/DigitBuilderConfig @@ -0,0 +1,11 @@ +# DigitBuilder config file + +verbosity 0 +ParametricModel 1 +#Reading in MC files +IsMC 1 +DigitChargeThr 20 +# There are three configurations: "PMT_only", "LAPPD_only", "All" +PhotoDetectorConfiguration All +#File must be in /pnfs/ space when loading on Fermilab cluster +LAPPDIDFile configfiles/EnergyReco/Train_Test/LAPPDIDs.txt diff --git a/configfiles/EnergyReco/Train_Test/EventSelectorConfig b/configfiles/EnergyReco/Train_Test/EventSelectorConfig new file mode 100644 index 000000000..ca6032811 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/EventSelectorConfig @@ -0,0 +1,15 @@ +# EventSelector config file + +verbosity 0 +MCPMTVolCut 0 +MCFVCut 1 +MCMRDCut 1 +MCPiKCut 0 +MCIsMuonCut 1 +MRDRecoCut 0 +RecoPMTVolCut 0 +RecoFVCut 0 +NHitCut 1 +PromptTrigOnly 1 +SaveStatusToStore 1 +Triggerword -1 diff --git a/configfiles/EnergyReco/Train_Test/FindMrdTracksConfig b/configfiles/EnergyReco/Train_Test/FindMrdTracksConfig new file mode 100644 index 000000000..7a9d11f7d --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/FindMrdTracksConfig @@ -0,0 +1,12 @@ +# FindMrdTracks Config File +# all variables retrieved with m_variables.Get() must be defined here! + +verbosity 0 +IsData 0 +OutputDirectory . +OutputFile STEC_MRDTracks_cluster40ns +DrawTruthTracks 0 # whether to add MC Truth track info for drawing in MrdPaddlePlot Tool + ## note you need to run that tool to actually view the tracks! +WriteTracksToFile 0 # should the track information be written to a ROOT-file? +SelectTriggerType 0 #should the loaded data be filtered by trigger type? +TriggerType Beam #options: Cosmic, Beam, No Loopback diff --git a/configfiles/EnergyReco/Train_Test/FindTrackLengthInWaterConfig b/configfiles/EnergyReco/Train_Test/FindTrackLengthInWaterConfig new file mode 100644 index 000000000..7744dae8c --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/FindTrackLengthInWaterConfig @@ -0,0 +1,6 @@ +#FindTrackLengthInWater Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbosity 1 +MaxTotalHitsToDNN 1100 +OutputDataFile data_for_trackLength_training.csv +DoTraining 1 diff --git a/configfiles/EnergyReco/Train_Test/HitCleanerConfig b/configfiles/EnergyReco/Train_Test/HitCleanerConfig new file mode 100644 index 000000000..36d7f0137 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/HitCleanerConfig @@ -0,0 +1,20 @@ +# HitCleaner config file + +verbosity 5 +Config 2 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +PmtMinPulseHeight 20 #minimum pulse height +PmtNeighbourRadius 60 #digit neighbouring distance [cm] +PmtMinNeighbourDigits 2 #minimum neighbour digits +PmtClusterRadius 60 #digit clustering distance [cm] +PmtTimeWindowN 10 #neighbouring time window [ns] +PmtTimeWindowC 10 #clustering time window [ns] +PmtMinHitsPerCluster 4 #number of hits per cluster +LappdMinPulseHeight 0 #minimum pulse height +LappdNeighbourRadius 25 #digit neighbouring distance [cm] +LappdMinNeighbourDigits 20 #minimum neighbour digits +LappdClusterRadius 25 #digit clustering distance [cm] +LappdTimeWindowN 1 #neighbouring time window [ns] +LappdTimeWindowC 1 #clustering time window [ns] +LappdMinHitsPerCluster 20 #number of digits per cluster +#MinClusterDigits 50 #minimum clustered digits (LAPPD+PMT) +MinClusterDigits 10 #minimum clustered digits (PMT-only) diff --git a/configfiles/EnergyReco/Train_Test/LAPPDIDs.txt b/configfiles/EnergyReco/Train_Test/LAPPDIDs.txt new file mode 100644 index 000000000..776b28912 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/LAPPDIDs.txt @@ -0,0 +1,5 @@ +11 +13 +14 +15 +17 diff --git a/configfiles/EnergyReco/Train_Test/LoadWCSimConfig b/configfiles/EnergyReco/Train_Test/LoadWCSimConfig new file mode 100644 index 000000000..001d2e929 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/LoadWCSimConfig @@ -0,0 +1,15 @@ +#LoadWCSim Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbose 0 + +InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standard/tank/pmt/wcsim_0* #James's files + +TriggerType Beam +SplitSubTriggers 1 + +WCSimVersion 10 ## should reflect the WCSim version of the files being loaded +HistoricTriggeroffset 0 ## time offset of digits relative to the trigger +#UseDigitSmearedTime 1 ## whether to use smeared digit time (T), or true time of first photon (F) +LappdNumStrips 60 ## num channels to construct from each LAPPD +LappdStripLength 100 ## relative x position of each LAPPD strip, for dual-sided readout [mm] +LappdStripSeparation 10 ## stripline separation, for calculating relative y position of each LAPPD strip [mm] diff --git a/configfiles/EnergyReco/Train_Test/LoadWCSimLAPPDConfig b/configfiles/EnergyReco/Train_Test/LoadWCSimLAPPDConfig new file mode 100644 index 000000000..2b82ff922 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/LoadWCSimLAPPDConfig @@ -0,0 +1,10 @@ +#LoadWCSimLAPPD Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbose 0 + +InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standard/tank/lappd/wcsim_lappd_0* #James's files + +WCSimVersion 10 ## should reflect the WCSim version of the files being loaded +InnerStructureRadius 1.3545 ## octagonal inner structure radius in m (from drawings 106.64") +HistoricTriggeroffset 0 +DrawDebugGraphs 0 ## whether to draw TPolyMarker3D's of hits diff --git a/configfiles/EnergyReco/Train_Test/MCParticlePropertiesConfig b/configfiles/EnergyReco/Train_Test/MCParticlePropertiesConfig new file mode 100644 index 000000000..752200268 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/MCParticlePropertiesConfig @@ -0,0 +1 @@ +verbosity 0 diff --git a/configfiles/EnergyReco/Train_Test/MCRecoEventLoaderConfig b/configfiles/EnergyReco/Train_Test/MCRecoEventLoaderConfig new file mode 100644 index 000000000..193bfd196 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/MCRecoEventLoaderConfig @@ -0,0 +1,12 @@ +# MCRecoEventLoader config file + +verbosity 0 +#Get Pion/Kaon counts from MC Truth +GetPionKaonInfo 1 +#Particle ID you want to load from parent particles; 13 for muon +ParticleID 13 +#coordinate shifts needed to put particles in tank origin coordinates +xshift 0.0 +yshift 14.46469 +zshift -168.1 + diff --git a/configfiles/EnergyReco/Train_Test/TimeClusteringConfig b/configfiles/EnergyReco/Train_Test/TimeClusteringConfig new file mode 100644 index 000000000..e3061f75f --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/TimeClusteringConfig @@ -0,0 +1,12 @@ +#TimeClustering config file + +verbosity 7 +MinDigitsForTrack 4 +MaxMrdSubEventDuration 40 +MinSubeventTimeSep 40 +MakeMrdDigitTimePlot 0 +LaunchTApplication 0 +IsData 0 +#OutputROOTFile TimeClustering_MRDTest28_cluster40ns +OutputROOTFile STEC_TimeClusteringOut +MapChankey_WCSimID ./configfiles/SimpleTankEnergyCalibrator/MRD_Chankey_WCSimID.dat diff --git a/configfiles/EnergyReco/Train_Test/ToolChainConfig b/configfiles/EnergyReco/Train_Test/ToolChainConfig new file mode 100644 index 000000000..59c288a77 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/ToolChainConfig @@ -0,0 +1,23 @@ +#ToolChain dynamic setup file + +##### Runtime Paramiters ##### +verbose 1 ## Verbosity level of ToolChain +error_level 0 # 0= do not exit, 1= exit on unhandeled errors only, 2= exit on unhandeled errors and handeled errors +attempt_recover 1 ## 1= will attempt to finalise if an execute fails + +###### Logging ##### +log_mode Interactive # Interactive=cout , Remote= remote logging system "serservice_name Remote_Logging" , Local = local file log; +log_local_path ./log +log_service LogStore + +###### Service discovery ##### Ignore these settings for local analysis +service_publish_sec -1 +service_kick_sec -1 + +##### Tools To Add ##### +Tools_File configfiles/EnergyReco/Train_Test/ToolsConfig ## list of tools to run and their config files + +##### Run Type ##### +Inline -1 ## number of Execute steps in program, -1 infinite loop that is ended by user +Interactive 0 ## set to 1 if you want to run the code interactively + diff --git a/configfiles/EnergyReco/Train_Test/ToolsConfig b/configfiles/EnergyReco/Train_Test/ToolsConfig new file mode 100644 index 000000000..8c68239ac --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/ToolsConfig @@ -0,0 +1,18 @@ +# train machine learning algorithm +myLoadWCSim LoadWCSim configfiles/EnergyReco/Train_Test/LoadWCSimConfig +myLoadWCSimLAPPD LoadWCSimLAPPD configfiles/EnergyReco/Train_Test/LoadWCSimLAPPDConfig +myMCParticleProperties MCParticleProperties configfiles/EnergyReco/Train_Test/MCParticlePropertiesConfig +MCRecoEventLoader MCRecoEventLoader configfiles/EnergyReco/Train_Test/MCRecoEventLoaderConfig +DigitBuilder DigitBuilder configfiles/EnergyReco/Train_Test/DigitBuilderConfig +HitCleaner HitCleaner configfiles/EnergyReco/Train_Test/HitCleanerConfig +ClusterFinder ClusterFinder configfiles/EnergyReco/Train_Test/ClusterFinderConfig +TimeClustering TimeClustering configfiles/EnergyReco/Train_Test/TimeClusteringConfig +EventSelector EventSelector configfiles/EnergyReco/Train_Test/EventSelectorConfig +FindMrdTracks FindMrdTracks configfiles/EnergyReco/Train_Test/FindMrdTracksConfig +VtxSeedGenerator VtxSeedGenerator configfiles/EnergyReco/Train_Test/VtxSeedGeneratorConfig +VtxSeedFineGrid VtxSeedFineGrid configfiles/EnergyReco/Train_Test/VtxSeedFineGridConfig +VtxExtendedVertexFinder VtxExtendedVertexFinder configfiles/EnergyReco/Train_Test/VtxExtendedVertexFinderConfig + +FindTrackLengthInWater FindTrackLengthInWater configfiles/EnergyReco/Train_Test/FindTrackLengthInWaterConfig +DNNTrackLengthTrainer PythonScript configfiles/EnergyReco/Train_Test/DNNTrackLengthTrainConfig +BDTMuonEnergyTrainer PythonScript configfiles/EnergyReco/Train_Test/BDTMuonEnergyTrainConfig diff --git a/configfiles/EnergyReco/Train_Test/VtxExtendedVertexFinderConfig b/configfiles/EnergyReco/Train_Test/VtxExtendedVertexFinderConfig new file mode 100644 index 000000000..59f74e5eb --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/VtxExtendedVertexFinderConfig @@ -0,0 +1,9 @@ +# VtxPointPositionFinder config file + +verbosity 5 + +# Only one of the options below should be set to 1. All the others should be set to 0 +UseTrueVertexAsSeed 0 +UsePointVertexAsSeed 0 +FitAllOnSeedGrid 1 +UseMeanTimeAsSeed 0 diff --git a/configfiles/EnergyReco/Train_Test/VtxSeedFineGridConfig b/configfiles/EnergyReco/Train_Test/VtxSeedFineGridConfig new file mode 100644 index 000000000..674e643a0 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/VtxSeedFineGridConfig @@ -0,0 +1,7 @@ +verbosity 5 +useDirectionGrid 0 +useSimpleDirection 1 +useTrueDir 0 +useMRDTrack 0 +usePastResolution 0 +multiGrid 1 diff --git a/configfiles/EnergyReco/Train_Test/VtxSeedGeneratorConfig b/configfiles/EnergyReco/Train_Test/VtxSeedGeneratorConfig new file mode 100644 index 000000000..014f6e4e7 --- /dev/null +++ b/configfiles/EnergyReco/Train_Test/VtxSeedGeneratorConfig @@ -0,0 +1,6 @@ +# VtxSeedGenerator config file + +verbosity 3 +NumberOfSeeds 300 +UseSeedGrid 1 +SeedType 2