Fixes to some scripts (#3)

* Add validate program for resolution

* Loop over all events to check momentum

* Update to simulation script

* Fix minor bug in script

* Possibly a quick generator for electron lund

* Adding different histograms

* Updates to sim and validate scripts

.gitignore

@@ -0,0 +1,7 @@
+
+*.ipynb
+*.o
+*.root
+release-validation-1/src/validate
+release-validation-1/test
+release-validation-1/validate

release-validation-1/gen.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+from __future__ import print_function
+from datetime import datetime
+import argparse
+import numpy as np
+import sys
+import os
+
+
+def gen_lund(args):
+    header = "           2  1.  1.  0 -1 0.0   0.0000  0.0000   0.0000   0.0000"
+    event1 = "    1  -1  2  11  0  0  0.0000  0.0000  {:.4f}  {:.4f}  0.0005  {:.4f}  {:.4f}  {:.4f}"
+    event2 = "    2   0  1  11  0  0  {:.4f}  {:.4f}  {:.4f}  {:.4f}  0.0005  {:.4f}  {:.4f}  {:.4f}"
+    energy = args.energy
+    theta = args.theta
+    phi = args.phi
+
+    with open(args.output_file, "w") as output:
+        for evt in range(0, args.num):
+            output.write(header)
+            output.write("\n")
+
+            if(args.sigma > 0):
+                energy = args.energy
+                theta = args.theta + \
+                    float(np.random.normal(0.0, args.sigma, 1))
+                phi = args.phi + float(np.random.normal(0.0, args.sigma, 1))
+
+            px = energy * np.sin(theta) * np.cos(phi)
+            py = energy * np.sin(theta) * np.sin(phi)
+            pz = energy * np.cos(theta)
+            vx = 0.0
+            vy = 0.0
+            vz = 0.0
+            E = np.sqrt(px * px + py * py + pz * pz)
+            output.write(event1.format(E, E, vx, vy, vz))
+            output.write("\n")
+            output.write(event2.format(px, py, pz, E, vx, vy, vz))
+            output.write("\n")
+        output.close()
+
+
+def main():
+    # Make argument parser
+    parser = argparse.ArgumentParser(
+        description="Basic generator for electrons")
+    parser.add_argument('-o', dest='output_file', type=str, nargs='?',
+                        help="Output filename", default='electron_gen.dat')
+    parser.add_argument('-E', dest='energy', type=float, nargs='?',
+                        help="Energy to generate electrons at", default=4.5)
+    parser.add_argument('-T', dest='theta', type=float, nargs='?',
+                        help="Theta angle to generate electrons at", default=20)
+    parser.add_argument('-P', dest='phi', type=float, nargs='?',
+                        help="phi angle to generate electrons at", default=20)
+    parser.add_argument('-s', dest='sigma', type=float, nargs='?',
+                        help="Generate electron with noise", default=0.0)
+    parser.add_argument('-N', dest='num', type=int, nargs='?',
+                        help="Number of electrons to generate", default=200)
+    args = parser.parse_args()
+
+    gen_lund(args)
+
+
+if __name__ == "__main__":
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("\n\nExiting")
+        sys.exit()

release-validation-1/run-clara-validate

@@ -0,0 +1,69 @@
+#!/usr/bin/env bash
+
+# Copyright (c) 2017.  Jefferson Lab (JLab). All rights reserved. Permission
+# to use, copy, modify, and distribute  this software and its documentation for
+# educational, research, and not-for-profit purposes, without fee and without a
+# signed licensing agreement.
+#
+# IN NO EVENT SHALL JLAB BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL
+# INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING
+# OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF JLAB HAS
+# BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+#
+# JLAB SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
+# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+# PURPOSE. THE CLARA SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY,
+# PROVIDED HEREUNDER IS PROVIDED "AS IS". JLAB HAS NO OBLIGATION TO PROVIDE
+# MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+#
+# This software was developed under the United States Government license.
+# For more information contact author at [email protected]
+# Department of Experimental Nuclear Physics, Jefferson Lab.
+
+BASEDIR=$(cd "$(dirname "${BASH_SOURCE[0]}")"/.. && pwd)
+
+if [ -z "${CLARA_HOME}" ]; then
+    CLARA_HOME=${BASEDIR}
+    export CLARA_HOME
+fi
+
+# find plugins directory
+if [ -n "${CLARA_PLUGINS}" ]; then
+    if [ ! -d "${CLARA_PLUGINS}" ]; then
+        echo "Error: \$CLARA_PLUGINS is not a directory."
+        exit 1
+    fi
+    plugins_dir="${CLARA_PLUGINS}"
+else
+    plugins_dir="${CLARA_HOME}/plugins"
+fi
+
+# set default classpath
+if [ -z "${CLASSPATH}" ]; then
+  CLASSPATH="${CLARA_HOME}/lib/*:${CLARA_HOME}/plugins/clas12/lib/services/*"
+  if [ -z ${CLAS12DIR+x} ]; then
+    for plugin in "${plugins_dir}"/*/; do
+      plugin=${plugin%*/}
+      if [ "${plugin##*/}" = "clas12" ]; then # COAT has special needs
+          CLASSPATH+=":${plugin}/lib/clas/*:${plugin}/lib/services/*"
+      else
+          CLASSPATH+=":${plugin}/services/*:${plugin}/lib/*"
+      fi
+    done
+  else
+      CLASSPATH+=":$CLAS12DIR/lib/clas/*:$CLAS12DIR/lib/services/*"
+  fi
+
+    CLASSPATH+=":${CLARA_HOME}/services/*"
+    export CLASSPATH
+fi
+
+if [ -n "${JAVA_OPTS}" ]; then
+    jvm_options=(${JAVA_OPTS})
+else
+    jvm_options=(-Xms256m)
+fi
+
+java_wrapper=${CLARA_HOME}/lib/clara/run-java
+
+exec "${java_wrapper}" "${jvm_options[@]}" org.jlab.clara.std.cli.ClaraShell "$@"

release-validation-1/sim.py

@@ -49,7 +49,7 @@ def split_lund(file_name, files):
     with open(file_name, "r") as text_file:
         lines = text_file.readlines()
         for line in lines:
-            if line.startswith("          "):
+            if line.startswith("        "):
                 num_events += 1
 
         out = 0
@@ -66,7 +66,7 @@ def split_lund(file_name, files):
                 else:
                     break
 
-            if line.startswith("          "):
+            if line.startswith("        "):
                 tot += 1
                 num_e = int(re.findall(r'\d+', line)[0]) + 1
                 for x in range(num_e):
@@ -75,21 +75,21 @@ def split_lund(file_name, files):
 
 
 def do_gemc(base):
-    cwd = os.getcwd()
+    base = os.getcwd() + "/" + base if base[0] != '/' else base
     with tempdir() as dir_temp:
         shutil.copy(base + ".dat", dir_temp + "/input.dat")
         with open(dir_temp + "/do_sim.sh", "w") as text_file:
-            text_file.write(r"""#!/bin/bash
-            source /jlab/2.2/ce/jlab.sh 2> /dev/null
-
-            /jlab/clas12Tags/4a.2.4/source/gemc /jlab/workdir/clas12.gcard -USE_GUI=0 -OUTPUT="evio, /jlab/workdir/shared/out.evio" -INPUT_GEN_FILE="LUND, /jlab/workdir/shared/input.dat"
-            """)
-
-        command = "docker run -v`pwd`:/jlab/workdir/shared --rm -it jeffersonlab/clas12tags:4a.2.4 bash /jlab/workdir/shared/do_sim.sh "
+            text_file.write("#!/bin/bash \n")
+            text_file.write("source /jlab/2.2/ce/jlab.sh \n")
+            text_file.write(
+                "/jlab/clas12Tags/4a.2.4/source/gemc /jlab/workdir/clas12.gcard -USE_GUI=0 -OUTPUT=\"evio, /jlab/workdir/shared/out.evio\" -INPUT_GEN_FILE=\"LUND, /jlab/workdir/shared/input.dat\" \n\n")
+            text_file.write("echo \"Gemc Exit Code: \" $? \n")
+            text_file.write("\n")
+
+        command = "docker run -v`pwd`:/jlab/workdir/shared --rm -it jeffersonlab/clas12tags:4a.2.4 bash /jlab/workdir/shared/do_sim.sh"
         out = " 1>" + base + ".out"
         err = " 2>" + base + ".err"
-        command = command + err + out
-        exit_code = os.system(command + " && exit 0")
+        exit_code = os.system(command + out + err)
         shutil.copy(dir_temp + "/out.evio", base + ".evio")
 
 
@@ -99,16 +99,25 @@ def main():
     parser.add_argument('-c', dest='cores', type=int, nargs='?',
                         help="Number of cores to use for simulation if not all the cores", default=0)
     parser.add_argument('-o', dest='output_dir', type=str, nargs='?',
-                        help="Output directory for final files", default=os.getcwd())
+                        help="Output directory for final files", default='.')
     parser.add_argument('-i', dest='events', type=str, nargs='?',
                         help="Input event file to run gemc over", default="11gev_sidis_500.dat")
+    parser.add_argument('-mc', dest='mc', action='store_true',
+                        help="Use gemc MC to generate electrons in a sector")
+    parser.set_defaults(mc=False)
+    parser.add_argument('-E', dest='energy', type=float, nargs='?',
+                        help="Energy for gemc MC to generate electrons at", default=4.5)
+    parser.add_argument('-T', dest='theta', type=float, nargs='?',
+                        help="Theta angle for gemc MC to generate electrons at", default=20)
+    parser.add_argument('-P', dest='phi', type=float, nargs='?',
+                        help="phi angle for gemc MC to generate electrons at", default=20)
 
     args = parser.parse_args()
 
-    if args.output_dir[-1] != '/':
-        args.output_dir = args.output_dir + '/'
-    if args.cores == 0 or args.cores > cpu_count():
-        args.cores = cpu_count()
+    args.output_dir = args.output_dir + \
+        '/' if args.output_dir[-1] != '/' else args.output_dir
+
+    args.cores = cpu_count() if args.cores == 0 or args.cores > cpu_count() else args.cores
 
     files = make_names(args.output_dir, args.events, args.cores)
     split_lund(args.events, files)

release-validation-1/src/Makefile

@@ -0,0 +1,21 @@
+ROOTLIBS	= $(shell root-config --libs)
+CXX = g++
+CXXFLAGS = -O3 -march=native -fPIC -w -g $(shell root-config --cflags)
+TARGET = validate
+
+MAINS = histogram.cpp main.cpp
+MAIN = $(MAINS:%.cpp=%.o)
+
+.PHONY: all clean
+
+all:	$(TARGET)
+
+$(MAIN): %.o : %.cpp
+		$(CXX) $(CXXFLAGS) -c $< -o $@
+
+$(TARGET): $(MAIN)
+	$(CXX) $(MAIN) $(CXXFLAGS) $(ROOTLIBS) -o $(TARGET)
+
+
+clean:
+	-rm -f $(TARGET) $(MAIN)

release-validation-1/src/colors.hpp

@@ -0,0 +1,22 @@
+#ifndef COLORS_H_GUARD
+#define COLORS_H_GUARD
+
+#define RESET "\033[0m"
+#define BLACK "\033[30m"              /* Black */
+#define RED "\033[31m"                /* Red */
+#define GREEN "\033[32m"              /* Green */
+#define YELLOW "\033[33m"             /* Yellow */
+#define BLUE "\033[34m"               /* Blue */
+#define MAGENTA "\033[35m"            /* Magenta */
+#define CYAN "\033[36m"               /* Cyan */
+#define WHITE "\033[37m"              /* White */
+#define BOLDBLACK "\033[1m\033[30m"   /* Bold Black */
+#define BOLDRED "\033[1m\033[31m"     /* Bold Red */
+#define BOLDGREEN "\033[1m\033[32m"   /* Bold Green */
+#define BOLDYELLOW "\033[1m\033[33m"  /* Bold Yellow */
+#define BOLDBLUE "\033[1m\033[34m"    /* Bold Blue */
+#define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */
+#define BOLDCYAN "\033[1m\033[36m"    /* Bold Cyan */
+#define BOLDWHITE "\033[1m\033[37m"   /* Bold White */
+
+#endif

release-validation-1/src/histogram.cpp

@@ -0,0 +1,66 @@
+/**************************************/
+/*																		*/
+/*  Created by Nick Tyler             */
+/*	University Of South Carolina      */
+/**************************************/
+#include "TCanvas.h"
+#include "histogram.hpp"
+
+Histogram::Histogram() {}
+
+Histogram::~Histogram() {}
+
+void Histogram::Fill_Res(double _px, double _py, double _pz, double _P, double _mc_px, double _mc_py, double _mc_pz,
+                         double _mc_P) {
+  momentum->Fill(_P);
+  momentum_x->Fill(_px);
+  momentum_y->Fill(_py);
+  momentum_z->Fill(_pz);
+
+  mom_rvg_x->Fill(_mc_px, _px);
+  mom_rvg_y->Fill(_mc_py, _py);
+  mom_rvg_z->Fill(_mc_pz, _pz);
+
+  resolution->Fill(_P - _mc_P);
+  resolution_x->Fill(_px - _mc_px);
+  resolution_y->Fill(_py - _mc_py);
+  resolution_z->Fill(_pz - _mc_pz);
+}
+
+void Histogram::Write() {
+  TCanvas c1("c1");
+  momentum->SetXTitle("momentum (GeV)");
+  momentum->Write();
+  momentum_x->SetXTitle("momentum (GeV)");
+  momentum_x->Write();
+  momentum_y->SetXTitle("momentum (GeV)");
+  momentum_y->Write();
+  momentum_z->SetXTitle("momentum (GeV)");
+  momentum_z->Write();
+
+  mom_rvg_x->SetXTitle("Gen X (GeV)");
+  mom_rvg_x->SetYTitle("Rec X (GeV)");
+  mom_rvg_x->SetOption("COLZ");
+  mom_rvg_x->Write();
+  mom_rvg_y->SetXTitle("Gen Y (GeV)");
+  mom_rvg_y->SetYTitle("Rec Y (GeV)");
+  mom_rvg_y->SetOption("COLZ");
+  mom_rvg_y->Write();
+  mom_rvg_z->SetXTitle("Gen Z (GeV)");
+  mom_rvg_z->SetYTitle("Rec Z (GeV)");
+  mom_rvg_z->SetOption("COLZ");
+  mom_rvg_z->Write();
+
+  resolution->SetXTitle("momentum (GeV)");
+  resolution->Fit("gaus", "QM+");
+  resolution->Write();
+  resolution_x->SetXTitle("momentum (GeV)");
+  resolution_x->Fit("gaus", "QM+");
+  resolution_x->Write();
+  resolution_y->SetXTitle("momentum (GeV)");
+  resolution_y->Fit("gaus", "QM+");
+  resolution_y->Write();
+  resolution_z->SetXTitle("momentum (GeV)");
+  resolution_z->Fit("gaus", "QM+");
+  resolution_z->Write();
+}

release-validation-1/src/histogram.hpp

@@ -0,0 +1,54 @@
+/************************************************************************/
+/*  Created by Nick Tyler*/
+/*	University Of South Carolina*/
+/************************************************************************/
+
+#ifndef HIST_H_GUARD
+#define HIST_H_GUARD
+#include "TF1.h"
+#include "TH1.h"
+#include "TH2.h"
+
+class Histogram {
+private:
+  int bins = 600;
+  double p_min = -6.0;
+  double p_max = 12.0;
+
+  TH1D *momentum =
+      new TH1D("ReconMom", "Reconstructed Momentum", bins, p_min, p_max);
+  TH1D *momentum_x =
+      new TH1D("ReconMom_x", "Reconstructed Momentum in x", bins, p_min, p_max);
+  TH1D *momentum_y =
+      new TH1D("ReconMom_y", "Reconstructed Momentum in y", bins, p_min, p_max);
+  TH1D *momentum_z =
+      new TH1D("ReconMom_z", "Reconstructed Momentum in z", bins, p_min, p_max);
+
+  TH2D *mom_rvg_x =
+      new TH2D("Recon_vs_gen_X", "Momentum X", bins, p_min,
+               p_max, bins, p_min, p_max);
+  TH2D *mom_rvg_y =
+      new TH2D("Recon_vs_gen_Y", "Momentum Y", bins, p_min,
+               p_max, bins, p_min, p_max);
+  TH2D *mom_rvg_z =
+      new TH2D("Recon_vs_gen_Z", "Momentum Z", bins, p_min,
+               p_max, bins, p_min, p_max);
+
+  TH1D *resolution = new TH1D("MomRes", "Momentum Resolution", 50, -2.5, 2.5);
+  TH1D *resolution_x =
+      new TH1D("MomRes_x", "Momentum Resolution in x", 50, -1.0, 1.0);
+  TH1D *resolution_y =
+      new TH1D("MomRes_y", "Momentum Resolution in y", 50, -1.0, 1.0);
+  TH1D *resolution_z =
+      new TH1D("MomRes_z", "Momentum Resolution in z", 50, -1.0, 1.0);
+
+public:
+  Histogram();
+  ~Histogram();
+
+  void Fill_Res(double px, double py, double pz, double P, double mc_px,
+                double mc_py, double mc_pz, double mc_P);
+  void Write();
+};
+
+#endif