diff --git a/benchmarks/roman_pots/Snakefile b/benchmarks/roman_pots/Snakefile new file mode 100644 index 00000000..e80af748 --- /dev/null +++ b/benchmarks/roman_pots/Snakefile @@ -0,0 +1,229 @@ +from itertools import product + +DETECTOR_PATH = os.environ["DETECTOR_PATH"] +DETECTOR_VERSION = re.search(r"epic-.*.-", DETECTOR_PATH).group().removeprefix("epic-").removesuffix("-") +SUBSYSTEM = "roman_pots" +BENCHMARK = "dense_neural_network" +DETECTOR_CONFIG = "epic_ip6" +NEVENTS_PER_FILE = 5 +NFILES = range(1,6) +MODEL_PZ = { + 'num_epochs' : [100], + 'learning_rate' : [0.01], + 'size_input' : [4], + 'size_output': [1], + 'n_layers' : [3,6], + 'size_first_hidden_layer' : [128], + 'multiplier' : [0.5], + 'leak_rate' : [0.025], +} +MODEL_PY = { + 'num_epochs' : [100], + 'learning_rate' : [0.01], + 'size_input' : [3], + 'size_output': [1], + 'n_layers' : [3,6], + 'size_first_hidden_layer' : [128], + 'multiplier' : [0.5], + 'leak_rate' : [0.025] +} +MODEL_PX = { + 'num_epochs' : [100], + 'learning_rate' : [0.01], + 'size_input' : [3], + 'size_output': [1], + 'n_layers' : [3,7], + 'size_first_hidden_layer' : [128], + 'multiplier' : [0.5], + 'leak_rate' : [0.025] +} + +rule all: + input: + expand("results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/raw_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_{index}.edm4hep.root", + index=NFILES), + expand("results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/processed_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_{index}.txt", + index=NFILES), + expand("results/"+DETECTOR_VERSION+"/{detector_config}/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_pz/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/model_pz.pt", + detector_config=DETECTOR_CONFIG, + num_epochs=MODEL_PZ["num_epochs"], + learning_rate=MODEL_PZ["learning_rate"], + size_input=MODEL_PZ["size_input"], + size_output=MODEL_PZ["size_output"], + n_layers=MODEL_PZ["n_layers"], + size_first_hidden_layer=MODEL_PZ["size_first_hidden_layer"], + multiplier=MODEL_PZ["multiplier"], + leak_rate=MODEL_PZ["leak_rate"] + ), + expand("results/"+DETECTOR_VERSION+"/{detector_config}/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_pz/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/LossVsEpoch_model_pz.png", + detector_config=DETECTOR_CONFIG, + num_epochs=MODEL_PZ["num_epochs"], + learning_rate=MODEL_PZ["learning_rate"], + size_input=MODEL_PZ["size_input"], + size_output=MODEL_PZ["size_output"], + n_layers=MODEL_PZ["n_layers"], + size_first_hidden_layer=MODEL_PZ["size_first_hidden_layer"], + multiplier=MODEL_PZ["multiplier"], + leak_rate=MODEL_PZ["leak_rate"] + ), + expand("results/"+DETECTOR_VERSION+"/{detector_config}/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_py/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/model_py.pt", + detector_config=DETECTOR_CONFIG, + num_epochs=MODEL_PY["num_epochs"], + learning_rate=MODEL_PY["learning_rate"], + size_input=MODEL_PY["size_input"], + size_output=MODEL_PY["size_output"], + n_layers=MODEL_PY["n_layers"], + size_first_hidden_layer=MODEL_PY["size_first_hidden_layer"], + multiplier=MODEL_PY["multiplier"], + leak_rate=MODEL_PY["leak_rate"] + ), + expand("results/"+DETECTOR_VERSION+"/{detector_config}/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_py/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/LossVsEpoch_model_py.png", + detector_config=DETECTOR_CONFIG, + num_epochs=MODEL_PY["num_epochs"], + learning_rate=MODEL_PY["learning_rate"], + size_input=MODEL_PY["size_input"], + size_output=MODEL_PY["size_output"], + n_layers=MODEL_PY["n_layers"], + size_first_hidden_layer=MODEL_PY["size_first_hidden_layer"], + multiplier=MODEL_PY["multiplier"], + leak_rate=MODEL_PY["leak_rate"] + ), + expand("results/"+DETECTOR_VERSION+"/{detector_config}/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_px/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/model_px.pt", + detector_config=DETECTOR_CONFIG, + num_epochs=MODEL_PX["num_epochs"], + learning_rate=MODEL_PX["learning_rate"], + size_input=MODEL_PX["size_input"], + size_output=MODEL_PX["size_output"], + n_layers=MODEL_PX["n_layers"], + size_first_hidden_layer=MODEL_PX["size_first_hidden_layer"], + multiplier=MODEL_PX["multiplier"], + leak_rate=MODEL_PX["leak_rate"] + ), + expand("results/"+DETECTOR_VERSION+"/{detector_config}/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_px/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/LossVsEpoch_model_px.png", + detector_config=DETECTOR_CONFIG, + num_epochs=MODEL_PX["num_epochs"], + learning_rate=MODEL_PX["learning_rate"], + size_input=MODEL_PX["size_input"], + size_output=MODEL_PX["size_output"], + n_layers=MODEL_PX["n_layers"], + size_first_hidden_layer=MODEL_PX["size_first_hidden_layer"], + multiplier=MODEL_PX["multiplier"], + leak_rate=MODEL_PX["leak_rate"] + ) + + + + +rule roman_pots_generate_events: + input: + script="steering_file.py" + params: + detector_path=DETECTOR_PATH, + nevents_per_file=NEVENTS_PER_FILE, + detector_config=DETECTOR_CONFIG + output: + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/raw_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_{index}.edm4hep.root" + shell: + """ + npsim --steeringFile {input.script} \ + --compactFile {params.detector_path}/{params.detector_config}.xml \ + --outputFile {output} \ + -N {params.nevents_per_file} + """ + +rule roman_pots_preprocess_model_training_data: + input: + data = "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/raw_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_{index}.edm4hep.root", + script = "preprocess_model_training_data.cxx" + output: + full = "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/processed_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_{index}.txt", + lo = "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/processed_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_lo_{index}.txt" + shell: + """ + root -q -b {input.script}\"(\\\"{input.data}\\\",\\\"{output.full}\\\",\\\"{output.lo}\\\")\" + """ + +rule roman_pots_train_model_pz: + input: + data = ["results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/processed_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_{index}.txt".format(index=index) for index in NFILES], + script = "train_dense_neural_network.py" + params: + detector_version=DETECTOR_VERSION, + detector_config=DETECTOR_CONFIG, + subsystem=SUBSYSTEM, + benchmark=BENCHMARK + output: + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_pz/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/model_pz.pt", + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_pz/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/LossVsEpoch_model_pz.png" + shell: + """ + python {input.script} --input_files {input.data} --model_name model_pz --model_dir results/{params.detector_version}/{params.detector_config}/detector_benchmarks/{params.subsystem}/{params.benchmark}/artifacts/model_pz/num_epochs_{wildcards.num_epochs}/learning_rate_{wildcards.learning_rate}/size_input_{wildcards.size_input}/size_output_{wildcards.size_output}/n_layers_{wildcards.n_layers}/size_first_hidden_layer_{wildcards.size_first_hidden_layer}/multiplier_{wildcards.multiplier}/leak_rate_{wildcards.leak_rate} --num_epochs {wildcards.num_epochs} --learning_rate {wildcards.learning_rate} --size_input {wildcards.size_input} --size_output {wildcards.size_output} --n_layers {wildcards.n_layers} --size_first_hidden_layer {wildcards.size_first_hidden_layer} --multiplier {wildcards.multiplier} --leak_rate {wildcards.leak_rate} + """ + +rule roman_pots_train_model_py: + input: + data = ["results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/processed_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_{index}.txt".format(index=index) for index in NFILES], + script = "train_dense_neural_network.py" + params: + detector_version=DETECTOR_VERSION, + detector_config=DETECTOR_CONFIG, + subsystem=SUBSYSTEM, + benchmark=BENCHMARK + output: + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_py/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/model_py.pt", + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_py/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/LossVsEpoch_model_py.png" + shell: + """ + python {input.script} --input_files {input.data} --model_name model_py --model_dir results/{params.detector_version}/{params.detector_config}/detector_benchmarks/{params.subsystem}/{params.benchmark}/artifacts/model_py/num_epochs_{wildcards.num_epochs}/learning_rate_{wildcards.learning_rate}/size_input_{wildcards.size_input}/size_output_{wildcards.size_output}/n_layers_{wildcards.n_layers}/size_first_hidden_layer_{wildcards.size_first_hidden_layer}/multiplier_{wildcards.multiplier}/leak_rate_{wildcards.leak_rate} --num_epochs {wildcards.num_epochs} --learning_rate {wildcards.learning_rate} --size_input {wildcards.size_input} --size_output {wildcards.size_output} --n_layers {wildcards.n_layers} --size_first_hidden_layer {wildcards.size_first_hidden_layer} --multiplier {wildcards.multiplier} --leak_rate {wildcards.leak_rate} + """ + +rule roman_pots_train_model_px: + input: + data = ["results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/processed_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_{index}.txt".format(index=index) for index in NFILES], + script = "train_dense_neural_network.py" + params: + detector_version=DETECTOR_VERSION, + detector_config=DETECTOR_CONFIG, + subsystem=SUBSYSTEM, + benchmark=BENCHMARK + output: + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_px/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/model_px.pt", + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_px/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/LossVsEpoch_model_px.png" + shell: + """ + python {input.script} --input_files {input.data} --model_name model_px --model_dir results/{params.detector_version}/{params.detector_config}/detector_benchmarks/{params.subsystem}/{params.benchmark}/artifacts/model_px/num_epochs_{wildcards.num_epochs}/learning_rate_{wildcards.learning_rate}/size_input_{wildcards.size_input}/size_output_{wildcards.size_output}/n_layers_{wildcards.n_layers}/size_first_hidden_layer_{wildcards.size_first_hidden_layer}/multiplier_{wildcards.multiplier}/leak_rate_{wildcards.leak_rate} --num_epochs {wildcards.num_epochs} --learning_rate {wildcards.learning_rate} --size_input {wildcards.size_input} --size_output {wildcards.size_output} --n_layers {wildcards.n_layers} --size_first_hidden_layer {wildcards.size_first_hidden_layer} --multiplier {wildcards.multiplier} --leak_rate {wildcards.leak_rate} + """ + +rule roman_pots_train_model_py_lo: + input: + data = ["results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/processed_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_lo_{index}.txt".format(index=index) for index in NFILES], + script = "train_dense_neural_network.py" + params: + detector_version=DETECTOR_VERSION, + detector_config=DETECTOR_CONFIG, + subsystem=SUBSYSTEM, + benchmark=BENCHMARK + output: + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_py/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/model_py_lo.pt", + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_py/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/LossVsEpoch_model_py_lo.png" + shell: + """ + python {input.script} --input_files {input.data} --model_name model_py --model_dir results/{params.detector_version}/{params.detector_config}/detector_benchmarks/{params.subsystem}/{params.benchmark}/artifacts/model_py/num_epochs_{wildcards.num_epochs}/learning_rate_{wildcards.learning_rate}/size_input_{wildcards.size_input}/size_output_{wildcards.size_output}/n_layers_{wildcards.n_layers}/size_first_hidden_layer_{wildcards.size_first_hidden_layer}/multiplier_{wildcards.multiplier}/leak_rate_{wildcards.leak_rate} --num_epochs {wildcards.num_epochs} --learning_rate {wildcards.learning_rate} --size_input {wildcards.size_input} --size_output {wildcards.size_output} --n_layers {wildcards.n_layers} --size_first_hidden_layer {wildcards.size_first_hidden_layer} --multiplier {wildcards.multiplier} --leak_rate {wildcards.leak_rate} + """ + +rule roman_pots_train_model_px_lo: + input: + data = ["results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/processed_data/"+DETECTOR_VERSION+"_"+DETECTOR_CONFIG+"_lo_{index}.txt".format(index=index) for index in NFILES], + script = "train_dense_neural_network.py" + params: + detector_version=DETECTOR_VERSION, + detector_config=DETECTOR_CONFIG, + subsystem=SUBSYSTEM, + benchmark=BENCHMARK + output: + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_px/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/model_px_lo.pt", + "results/"+DETECTOR_VERSION+"/"+DETECTOR_CONFIG+"/detector_benchmarks/"+SUBSYSTEM+"/"+BENCHMARK+"/artifacts/model_px/num_epochs_{num_epochs}/learning_rate_{learning_rate}/size_input_{size_input}/size_output_{size_output}/n_layers_{n_layers}/size_first_hidden_layer_{size_first_hidden_layer}/multiplier_{multiplier}/leak_rate_{leak_rate}/LossVsEpoch_model_px_lo.png" + shell: + """ + python {input.script} --input_files {input.data} --model_name model_px --model_dir results/{params.detector_version}/{params.detector_config}/detector_benchmarks/{params.subsystem}/{params.benchmark}/artifacts/model_px/num_epochs_{wildcards.num_epochs}/learning_rate_{wildcards.learning_rate}/size_input_{wildcards.size_input}/size_output_{wildcards.size_output}/n_layers_{wildcards.n_layers}/size_first_hidden_layer_{wildcards.size_first_hidden_layer}/multiplier_{wildcards.multiplier}/leak_rate_{wildcards.leak_rate} --num_epochs {wildcards.num_epochs} --learning_rate {wildcards.learning_rate} --size_input {wildcards.size_input} --size_output {wildcards.size_output} --n_layers {wildcards.n_layers} --size_first_hidden_layer {wildcards.size_first_hidden_layer} --multiplier {wildcards.multiplier} --leak_rate {wildcards.leak_rate} + """ diff --git a/benchmarks/roman_pots/preprocess_model_training_data.cxx b/benchmarks/roman_pots/preprocess_model_training_data.cxx new file mode 100644 index 00000000..5d327ac7 --- /dev/null +++ b/benchmarks/roman_pots/preprocess_model_training_data.cxx @@ -0,0 +1,161 @@ +//------------------------- +// +// Hit reader to relate hits at Roman Pots to momentum vectors from MC. +// +// Input(s): output file from npsim particle gun for RP particles. +// +// Output(s): txt file with training information with px_mc, py_mc, pz_mc, x_rp, slope_xrp, y_rp, slope_yrp +// +// +// Author: Alex Jentsch +//------------------------ +//Low PT preprocessing added by David Ruth + +using namespace std; + +void preprocess_model_training_data(TString inputFile, TString outputFile, TString outputFile_lo) { + + string fileName; + TFile* inputRootFile; + TTree* rootTree; + + // MC information + TH1D* h_eta_MC = new TH1D("h_eta", ";Pseudorapidity, #eta", 100, 0.0, 15.0); + TH1D* h_px_MC = new TH1D("px_MC", ";p_{x} [GeV/c]", 100, -10.0, 10.0); + TH1D* h_py_MC = new TH1D("py_MC", ";p_{y} [GeV/c]", 100, -10.0, 10.0); + TH1D* h_pt_MC = new TH1D("pt_MC", ";p_{t} [GeV/c]", 100, 0.0, 2.0); + TH1D* h_pz_MC = new TH1D("pz_MC", ";p_{z} [GeV/c]", 100, 0.0, 320.0); + TH1D* h_theta_MC = new TH1D("theta_MC", ";#theta [mrad]", 100, 0.0, 25.0); + TH1D* h_phi_MC = new TH1D("phi_MC", ";#phi [rad]", 100, -3.2, 3.2); + + // Roman pots + TH1D* h_px_RomanPots = new TH1D("px_RomanPots", ";p_{x} [GeV/c]", 100, -10.0, 10.0); + TH1D* h_py_RomanPots = new TH1D("py_RomanPots", ";p_{y} [GeV/c]", 100, -10.0, 10.0); + TH1D* h_pt_RomanPots = new TH1D("pt_RomanPots", ";p_{t} [GeV/c]", 100, 0.0, 2.0); + TH1D* h_pz_RomanPots = new TH1D("pz_RomanPots", ";p_{z} [GeV/c]", 100, 0.0, 320.0); + + int fileCounter = 0; + int iEvent = 0; + + inputRootFile = new TFile(inputFile); + if (!inputRootFile) { + cout << "MISSING_ROOT_FILE" << fileName << endl; + return; + } + + TTree* evtTree = (TTree*)inputRootFile->Get("events"); + + int numEvents = evtTree->GetEntries(); + + TTreeReader tree_reader(evtTree); // !the tree reader + + ofstream outputTrainingFile; + outputTrainingFile.open(outputFile); + + ofstream outputTrainingFile_lo; + outputTrainingFile_lo.open(outputFile_lo); + // MC particles + TTreeReaderArray mc_px_array = { tree_reader, "MCParticles.momentum.x" }; + TTreeReaderArray mc_py_array = { tree_reader, "MCParticles.momentum.y" }; + TTreeReaderArray mc_pz_array = { tree_reader, "MCParticles.momentum.z" }; + TTreeReaderArray mc_mass_array = { tree_reader, "MCParticles.mass" }; + TTreeReaderArray mc_pdg_array = { tree_reader, "MCParticles.PDG" }; + TTreeReaderArray mc_genStatus_array = { tree_reader, "MCParticles.generatorStatus" }; + + // Roman pots -- momentum vector + + // hit locations (for debugging) + TTreeReaderArray global_hit_RP_x = { tree_reader, "ForwardRomanPotHits.position.x" }; + TTreeReaderArray global_hit_RP_y = { tree_reader, "ForwardRomanPotHits.position.y" }; + TTreeReaderArray global_hit_RP_z = { tree_reader, "ForwardRomanPotHits.position.z" }; + + cout << "file has " << evtTree->GetEntries() << " events..." << endl; + + tree_reader.SetEntriesRange(0, evtTree->GetEntries()); + while (tree_reader.Next()) { + + cout << "Reading event: " << iEvent << endl; + + if (mc_px_array.GetSize() > 1) { + cout << "Event has more than one particle -- skip..." << endl; + } + + bool hasMCProton = false; + bool hitLayerOne = false; + bool hitLayerTwo = false; + bool lowPT = false; + + double mcProtonMomentum[3]; + TVector3 mctrk; + + for (int imc = 0; imc < mc_px_array.GetSize(); imc++) { + mctrk.SetXYZ(mc_px_array[imc], mc_py_array[imc], mc_pz_array[imc]); + + if (mc_pdg_array[imc] == 2212 && mc_genStatus_array[imc] == 1) { //only checking for protons here -- change as desired + + mcProtonMomentum[0] = mctrk.Px(); + mcProtonMomentum[1] = mctrk.Py(); + mcProtonMomentum[2] = mctrk.Pz(); + hasMCProton = true; + } + } + + double hit1minZ = 32547.3-923/2.0; + double hit1maxZ = 32547.3+923/2.0; + double hit2minZ = 34245.5-923/2.0; + double hit2maxZ = 34245.5+923/2.0; + + double rpHitLayerOne[3]; + double rpHitLayerTwo[3]; + + double slopeXRP = 0.0; + double slopeYRP = 0.0; + + // roman pots reco tracks + for (int iRPPart = 0; iRPPart < global_hit_RP_x.GetSize(); iRPPart++) { + + if (global_hit_RP_z[iRPPart] > hit1minZ && global_hit_RP_z[iRPPart] < hit1maxZ) { + + rpHitLayerOne[0] = global_hit_RP_x[iRPPart]; + rpHitLayerOne[1] = global_hit_RP_y[iRPPart]; + rpHitLayerOne[2] = global_hit_RP_z[iRPPart]; + hitLayerOne = true; + } + + if (global_hit_RP_z[iRPPart] > hit2minZ && global_hit_RP_z[iRPPart] < hit2maxZ) { + + rpHitLayerTwo[0] = global_hit_RP_x[iRPPart]; + rpHitLayerTwo[1] = global_hit_RP_y[iRPPart]; + rpHitLayerTwo[2] = global_hit_RP_z[iRPPart]; + hitLayerTwo = true; + } + } + if (hasMCProton) { + double Pt = sqrt(pow(mcProtonMomentum[0],2) + pow(mcProtonMomentum[1],2)); + if(Pt < 0.3) { + lowPT = true; + } + } + if (hasMCProton && hitLayerOne && hitLayerTwo) { + + double slope_x = (rpHitLayerTwo[0] - rpHitLayerOne[0]) / (rpHitLayerTwo[2] - rpHitLayerOne[2]); + double slope_y = (rpHitLayerTwo[1] - rpHitLayerOne[1]) / (rpHitLayerTwo[2] - rpHitLayerOne[2]); + + outputTrainingFile << mcProtonMomentum[0] << "\t" << mcProtonMomentum[1] << "\t" << mcProtonMomentum[2] << "\t"; + outputTrainingFile << rpHitLayerTwo[0] << "\t" << slope_x << "\t" << rpHitLayerTwo[1] << "\t" << slope_y << endl; + if (lowPT) { + outputTrainingFile_lo << mcProtonMomentum[0] << "\t" << mcProtonMomentum[1] << "\t" << mcProtonMomentum[2] << "\t"; + outputTrainingFile_lo << rpHitLayerTwo[0] << "\t" << slope_x << "\t" << rpHitLayerTwo[1] << "\t" << slope_y << endl; + } + } + + iEvent++; + } // event loop + + inputRootFile->Close(); + outputTrainingFile.close(); + outputTrainingFile_lo.close(); + + return; +} + diff --git a/benchmarks/roman_pots/steering_file.py b/benchmarks/roman_pots/steering_file.py new file mode 100644 index 00000000..2ec59b01 --- /dev/null +++ b/benchmarks/roman_pots/steering_file.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python +""" +DD4hep simulation with some argument parsing +Based on M. Frank and F. Gaede runSim.py + @author A.Sailer + @version 0.1 +Modified with standard EIC EPIC requirements. +""" +from __future__ import absolute_import, unicode_literals +import logging +import sys +import math + +#from DDSim import SystemOfUnits +from DDSim.DD4hepSimulation import DD4hepSimulation +from g4units import mm, GeV, MeV, radian + +SIM = DD4hepSimulation() +################################################################ +# Steering file to shoot particles into the Roman Pots +# used to generate the training data for RP ML algorithm. +# +# The default magnetic field settings are for 18x275 GeV. +################################################################# + +############################################################# +# Particle Gun simulations use these options +# --> comment out if using MC input +############################################################# + +#SIM.enableDetailedShowerMode = False +SIM.enableG4GPS = False +SIM.enableG4Gun = False +SIM.enableGun = True +#SIM.gun.direction = (math.sin(-0.025*radian), 0, math.cos(-0.025*radian)) +SIM.crossingAngleBoost = -0.025*radian + + +## InputFiles for simulation .stdhep, .slcio, .HEPEvt, .hepevt, .hepmc, .pairs files are supported +SIM.inputFiles = [] +## Macro file to execute for runType 'run' or 'vis' +SIM.macroFile = "" +## number of events to simulate, used in batch mode +##SIM.numberOfEvents = 5000 +## Outputfile from the simulation,only lcio output is supported +##SIM.outputFile = "RP_eDep_and_thresholds_10_6_2023_1k_events.edm4hep.root" +SIM.runType = "run" +## FourVector of translation for the Smearing of the Vertex position: x y z t +SIM.vertexOffset = [0.0, 0.0, 0.0, 0.0] +## FourVector of the Sigma for the Smearing of the Vertex position: x y z t +SIM.vertexSigma = [0.0, 0.0, 0.0, 0.0] + + + +################################################################################ +## Configuration for the magnetic field (stepper) +################################################################################ +SIM.field.delta_chord = 1e-03 +SIM.field.delta_intersection = 1e-03 +SIM.field.delta_one_step = .5e-01*mm +SIM.field.eps_max = 1e-03 +SIM.field.eps_min = 1e-04 +#SIM.field.equation = "Mag_UsualEqRhs" +SIM.field.largest_step = 100.0*mm +SIM.field.min_chord_step = 1.e-2*mm +SIM.field.stepper = "HelixSimpleRunge" + +## choose the distribution of the random direction for theta +## +## Options for random distributions: +## +## 'uniform' is the default distribution, flat in theta +## 'cos(theta)' is flat in cos(theta) +## 'eta', or 'pseudorapidity' is flat in pseudorapity +## 'ffbar' is distributed according to 1+cos^2(theta) +## +## Setting a distribution will set isotrop = True +## +SIM.gun.distribution = 'uniform' +SIM.gun.energy = 275.0*GeV ## default energy value is MeV +#SIM.gun.momentumMin = 274.0*GeV +#SIM.gun.momentumMax = 275.0*GeV + +## isotropic distribution for the particle gun +## +## use the options phiMin, phiMax, thetaMin, and thetaMax to limit the range of randomly distributed directions +## if one of these options is not None the random distribution will be set to True and cannot be turned off! +## +SIM.gun.isotrop = True +#SIM.gun.multiplicity = 1 +SIM.gun.particle = "proton" + +## Minimal azimuthal angle for random distribution +##SIM.gun.phiMin = 0 + +## position of the particle gun, 3 vector. unit mm +##SIM.gun.position = (0, 0, 0) +SIM.gun.thetaMin = 0.000*radian +SIM.gun.thetaMax = 0.005*radian + + +################################################################################ +## Configuration for the output levels of DDG4 components +################################################################################ + +## Output level for input sources +SIM.output.inputStage = 3 + +## Output level for Geant4 kernel +SIM.output.kernel = 3 + +## Output level for ParticleHandler +SIM.output.part = 3 + +## Output level for Random Number Generator setup +SIM.output.random = 6 + + +################################################################################ +## Configuration for the Particle Handler/ MCTruth treatment +################################################################################ + +## Enable lots of printout on simulated hits and MC-truth information +SIM.part.enableDetailedHitsAndParticleInfo = False + +## Keep all created particles +SIM.part.keepAllParticles = True + +## Minimal distance between particle vertex and endpoint of parent after +## which the vertexIsNotEndpointOfParent flag is set +## +SIM.part.minDistToParentVertex = 2.2e-14 + +## MinimalKineticEnergy to store particles created in the tracking region +#SIM.part.minimalKineticEnergy = 1*MeV + +## Printout at End of Tracking +SIM.part.printEndTracking = False + +## Printout at Start of Tracking +SIM.part.printStartTracking = False + +## List of processes to save, on command line give as whitespace separated string in quotation marks +SIM.part.saveProcesses = ['Decay'] + + +################################################################################ +## Configuration for the PhysicsList +################################################################################ +SIM.physics.decays = True +SIM.physics.list = "FTFP_BERT" # "FTFP_BERT" + +################################################################################ +## Properties for the random number generator +################################################################################ + +## If True, calculate random seed for each event based on eventID and runID +## allows reproducibility even when SkippingEvents +SIM.random.enableEventSeed = False +SIM.random.file = None +SIM.random.luxury = 1 +SIM.random.replace_gRandom = True +SIM.random.seed = None +SIM.random.type = None diff --git a/benchmarks/roman_pots/train_dense_neural_network.py b/benchmarks/roman_pots/train_dense_neural_network.py new file mode 100644 index 00000000..36bd9bd6 --- /dev/null +++ b/benchmarks/roman_pots/train_dense_neural_network.py @@ -0,0 +1,163 @@ +import pandas as pd +import numpy as np +import seaborn as sb +import torch +import torch.nn as nn +import torch.optim as optim +import torch.optim.lr_scheduler as lr_scheduler +import matplotlib.pyplot as plt +import argparse +import sys +import hashlib + +torch.set_default_dtype(torch.float32) + +if torch.cuda.is_available(): + device = torch.device("cuda") + print("GPU is available!") +else: + device = torch.device("cpu") + print("GPU not found. Using CPU.") + +class NeuralNet(nn.Module): + def __init__(self, size_input, size_output, n_layers, size_first_hidden_layer=128, multiplier=0.5, leak_rate=0.025): + super().__init__() + self.layers = nn.ModuleList() + + size_current_hidden_layer = size_first_hidden_layer + + self.layers.append(nn.Linear(size_input, size_current_hidden_layer)) + for i in range(n_layers - 2): + self.layers.append(nn.LeakyReLU(leak_rate)) + self.layers.append(nn.Linear(size_current_hidden_layer, int(size_current_hidden_layer * multiplier))) + size_current_hidden_layer = int(size_current_hidden_layer * multiplier) + self.layers.append(nn.LeakyReLU(leak_rate)) + self.layers.append(nn.Linear(size_current_hidden_layer, size_output)) + + print("Create a network with the following layers:") + for layer in self.layers: + print(layer) + + def forward(self, x): + for layer in self.layers: + x = layer(x) + return x + +def standardize(x): + mean = torch.mean(x, axis=0) + std = torch.std(x, axis=0) + standardized_tensor = (x - mean) / std + return standardized_tensor, mean, std + +def train_model(input_tensor, target_tensor, model, hyperparameters): + # Send model to device + model=model.to(device) + + # Define the loss function and optimizer + criterion = torch.nn.HuberLoss(reduction='mean', delta=1.0) + optimizer = torch.optim.Adam(model.parameters(), lr=hyperparameters.learning_rate) + + # Create a learning rate scheduler + scheduler = lr_scheduler.ReduceLROnPlateau(optimizer,'min',patience=100,cooldown=100,factor=0.5,threshold=1e-4,verbose=True) + + # Track the losses + losses = [] + + # Train the model + for epoch in range(hyperparameters.num_epochs): + # Forward pass + inputs, targets = input_tensor.to(device), target_tensor.to(device) + predictions = model(inputs) + loss = criterion(predictions, targets) + + # Backward and optimize + optimizer.zero_grad() + loss.backward() + optimizer.step() + + # Track the loss value at each epoch + losses.append(loss.item()) + + # Step the learning rate scheduler + scheduler.step(loss) + + # Print progress + if (epoch + 1) % 10 == 0: + print("Epoch "+str(epoch+1)+"/"+str(hyperparameters.num_epochs)+", Loss: "+"{0:0.10f}".format(loss.item())) + + # Plot the loss values + plt.figure() + plt.plot(range(1, hyperparameters.num_epochs+1), losses) + plt.xlabel('Epoch') + plt.ylabel('Loss') + plt.title('Loss as a Function of Epoch') + plt.yscale('log') + plt.savefig(hyperparameters.model_dir+"/LossVsEpoch_"+hyperparameters.model_name+".png") + + torch.jit.script(model).save(hyperparameters.model_dir+"/"+hyperparameters.model_name+".pt") + return + +def run_experiment(hyperparameters): + + # Load training data in tensors + training_data = pd.DataFrame() + + for i in hyperparameters.input_files: + temp_training_data = pd.read_csv(i, delimiter='\t', header=None) + training_data = pd.concat([training_data, temp_training_data], ignore_index=True) + + training_RP_pos_tensor = torch.tensor(training_data.iloc[:,3:7].values, dtype=torch.float32) + training_MC_mom_tensor = torch.tensor(training_data.iloc[:,0:3].values, dtype=torch.float32) + + # Standardize training data + match hyperparameters.model_name: + case "model_pz": + source = training_RP_pos_tensor + scaled_source, mean_source, std_source = standardize(source) + target = training_MC_mom_tensor[:,2].unsqueeze(1) + + case "model_py": + source = torch.cat((training_RP_pos_tensor[:,2:4], training_MC_mom_tensor[:,2].unsqueeze(1)), 1) + scaled_source, mean_source, std_source = standardize(source) + target = training_MC_mom_tensor[:,1].unsqueeze(1) + + case "model_px": + source = torch.cat((training_RP_pos_tensor[:,0:2], training_MC_mom_tensor[:,2].unsqueeze(1)), 1) + scaled_source, mean_source, std_source = standardize(source) + target = training_MC_mom_tensor[:,0].unsqueeze(1) + + case _: + print("No model name provided. Stop further processing") + return + + # Initialize models + initial_model = NeuralNet(size_input=int(hyperparameters.size_input), + size_output=int(hyperparameters.size_output), + n_layers=int(hyperparameters.n_layers), + size_first_hidden_layer=int(hyperparameters.size_first_hidden_layer), + multiplier=float(hyperparameters.multiplier), + leak_rate=float(hyperparameters.leak_rate)) + + # Train models + train_model(scaled_source, target, initial_model, hyperparameters) + + # Print end statement + print("Training completed using "+str(len(hyperparameters.input_files))+" files with "+str(training_RP_pos_tensor.shape[0])+" eligible events") + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Train neural network model for roman pots") + parser.add_argument('--input_files', type=str, nargs='+', required=True, help='Specify a location of input files.') + parser.add_argument('--model_name', type=str, required=True, help='Specify model name.') + parser.add_argument('--model_dir', type=str, required=True, help='Specify location to save model') + parser.add_argument('--num_epochs', type=int, required=True, help='Specify number of epochs') + parser.add_argument('--learning_rate', type=float, required=True, help='Specify learning rate') + parser.add_argument('--size_input', type=int, required=True, help='Specify input size') + parser.add_argument('--size_output', type=int, required=True, help='Specify output size') + parser.add_argument('--n_layers', type=int, required=True, help='Specify number of layers') + parser.add_argument('--size_first_hidden_layer', type=int, required=True, help='Size of first hidden layer') + parser.add_argument('--multiplier', type=float, required=True, help='Specify mutilplier to calculate size of subsequent hidden layers') + parser.add_argument('--leak_rate', type=float, required=True, help='Specify leak rate') + hyperparameters = parser.parse_args() + run_experiment(hyperparameters) + +