long time, no commit

Author: RaulRPrado

changeGroup

@@ -0,0 +1,50 @@
+#!/bin/bash
+
+set -o errexit   # abort on nonzero exitstatus
+set -o pipefail  # don't hide errors within pipes
+
+echo "Starting changeGroup"
+
+# shellcheck disable=SC1091,SC1090
+source "util.sh"
+
+atm="summer"
+zenith=""
+wobble=""
+n_min=0        # min run
+n_max=200001     # max run (100000)
+source_dir=${PWD} 
+readonly runs_per_job=100
+readonly do_submit="true"
+readonly force="false"
+readonly n_mail=5
+
+# shellcheck disable=SC2068
+collect_arguments 4 atm zenith wobble n_max $@
+
+all_n=$(seq "${n_min}" "${n_max}")
+
+atm=$(validate_atm "${atm}")
+zenith=$(validate_zenith "${zenith}")
+all_wobble="0.0 0.25 0.5 0.75 1.0 1.25 1.5 1.75 2.0 2.8"
+
+# shellcheck disable=SC2086
+for wobble in ${all_wobble}; do 
+    echo "Starting wobble ${wobble}"
+
+    gro_file=$(groptics_file "10" "${zenith}" "${atm}" "${wobble}")
+    gro_dir=$(dirname "${gro_file}")
+    if [ ! -d "${gro_dir}" ]; then
+        echo "Directory ${gro_dir} does not exist - skipping"
+        continue
+    fi
+
+    for n in ${all_n}; do
+        run=$(compute_run "${zenith}" "${n}")
+        gro_file=$(groptics_file "${run}" "${zenith}" "${atm}" "${wobble}")
+        if [ -f "${gro_file}" ]; then
+            chgrp veritas "${gro_file}"
+        fi
+    done
+
+done

checkShowerPars.py

@@ -0,0 +1,114 @@
+#!/usr/bin/python3
+
+import matplotlib.pyplot as plt
+import numpy as np
+import argparse
+import math
+import os.path
+from matplotlib.backends.backend_pdf import PdfPages
+
+import ROOT
+
+labelsize = 15
+plt.rc('font', family='serif', size=labelsize)
+plt.rc('xtick', labelsize=labelsize)
+plt.rc('ytick', labelsize=labelsize)
+plt.rc('text', usetex=True)
+
+
+def check_woff(atm, zenith, wobble, nsb, pdf):
+    print('check_woff({}, {}, {}, {})'.format(atm, zenith, wobble, nsb))
+
+    woff_rad = list()
+    for i_split in range(10):
+        run = str(961201 + i_split)
+        try:
+            filename = (
+                '/lustre/fs23/group/veritas/V6_DESY/EffAreaTesting/v501/CARE/V6_ATM'
+                + str(atm) + '_gamma_ze' + str(zenith)
+                + '_TL5035MA20/ze' + str(zenith) + 'deg_offset' + str(wobble) + 'deg_NSB'
+                + str(nsb) + 'MHz/' + run + '.root'
+            )
+
+            if not os.path.isfile(filename):
+                # print('File {} does not exists'.format(filename))
+                continue
+
+            try:
+                file = ROOT.TFile.Open(filename)
+            except:
+                print('Could not open file {}'.format(filename))
+                continue
+
+            for entry in file.showerpars:
+                if float(entry.Xoff[0]) > -99 and not float(entry.Xoff[0]) == 0.:
+                    woff_rad.append(math.sqrt(entry.Xoff[0]**2 + entry.Yoff[0]**2))
+        except:
+            print('Could not read file for run {}'.format(run))
+            continue
+
+    mean_woff = np.sum(woff_rad) / len(woff_rad)
+    diff = math.fabs(mean_woff - float(wobble))
+    sigma = np.std(woff_rad)
+
+    fig = plt.figure(figsize=(8, 6), tight_layout=True)
+
+    ax = plt.gca()
+    ax.set_title('NSB = {}, mu = {:.3f}, sigma = {:.3f}'.format(nsb, mean_woff, sigma))
+    ax.set_xlabel('Woff')
+
+    ax.hist(woff_rad, bins=np.linspace(0.0, 2.5, 100))
+
+    pdf.savefig(fig)
+    plt.close()
+
+    return
+
+
+if __name__ == '__main__':
+
+    parser = argparse.ArgumentParser(
+        description='Plot Eff Area for a certain zenith angle simulation set.'
+    )
+    parser.add_argument(
+        '-z',
+        '--zenith',
+        help='Zenith angle (in deg)',
+        type=str,
+        required=True,
+        choices=['00', '20', '30', '35', '40', '45', '50', '55', '60']
+    )
+    parser.add_argument(
+        '-a',
+        '--atm',
+        help='Atm (61 or 62)',
+        type=str,
+        required=True,
+        choices=['61', '62']
+    )
+    parser.add_argument(
+        '-w',
+        '--wobble',
+        help='Wobble offset',
+        type=float,
+        required=False,
+        choices=[0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0]
+    )
+    zenith = parser.parse_args().zenith
+    atm = parser.parse_args().atm
+    wob = parser.parse_args().wobble
+
+    if wob is None:
+        all_wobble = [0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0]
+        figName = 'figures/Woff_atm{}_ze{}.pdf'.format(atm, zenith)
+    else:
+        all_wobble = [wob]
+        figName = 'figures/Woff_atm{}_ze{}_wob{}.pdf'.format(atm, zenith, wob)
+
+    all_nsb = [50, 75, 100, 130, 160, 200, 250, 300, 350, 400, 450]
+    # all_nsb = [400]
+
+    with PdfPages(figName) as pdf:
+        for nsb in all_nsb:
+            for wob in all_wobble:
+                check_woff(atm, zenith, wob, nsb, pdf)

compareEffArea.py

@@ -0,0 +1,140 @@
+#!/usr/bin/python3
+
+import matplotlib.pyplot as plt
+import numpy as np
+import argparse
+import os.path
+import sys
+import math
+from matplotlib.backends.backend_pdf import PdfPages
+
+import ROOT
+
+labelsize = 15
+plt.rc('font', family='serif', size=labelsize)
+plt.rc('xtick', labelsize=labelsize)
+plt.rc('ytick', labelsize=labelsize)
+plt.rc('text', usetex=True)
+
+
+def get_eff_area(filename, azimuth=0, index=2.0):
+    lge, eff, eff_err = list(), list(), list()
+
+    if not os.path.isfile(filename):
+        print('File {} does not exists'.format(filename))
+        return None, None, None
+
+    file = ROOT.TFile.Open(filename)
+    for entry in file.fEffArea:
+        if (entry.az == azimuth and entry.index == index):
+            lge = list(entry.e0)
+            eff = list(entry.eff)
+            eff_err = list(entry.eff_error)
+            return lge, eff, eff_err
+        else:
+            continue
+    print('get_eff_area failed')
+    return None, None, None
+
+
+if __name__ == '__main__':
+
+    if len(sys.argv) != 3:
+        print('Wrong number of args')
+        exit()
+
+    lge0, eff0, eff_err0 = get_eff_area(sys.argv[1])
+    lge1, eff1, eff_err1 = get_eff_area(sys.argv[2])
+
+    print(lge0)
+    print(lge1)
+
+    diff = list()
+    rel_err0, rel_err1 = list(), list()
+    for i in range(len(lge0)):
+        d = (eff0[i] - eff1[i]) / math.sqrt(eff_err0[i]*eff_err0[i] + eff_err1[i]*eff_err1[i])
+        diff.append(d)
+        rel_err0.append(eff_err0[i] / eff0[i])
+        rel_err1.append(eff_err1[i] / eff1[i])
+
+    figName = 'figures/EffArea_testingSplit.pdf'
+    with PdfPages(figName) as pdf:
+        fig = plt.figure(figsize=(8, 6), tight_layout=True)
+        ax = plt.gca()
+        ax.set_yscale('log')
+        ax.set_xlabel(r'log$_{10}$($E$/TeV)')
+        ax.set_ylabel(r'$A_\mathrm{eff}$ (cm$^2$)')
+
+        ax.errorbar(
+            lge0,
+            eff0,
+            yerr=eff_err0,
+            marker='o',
+            linestyle='none'
+        )
+        ax.errorbar(
+            lge1,
+            eff1,
+            yerr=eff_err1,
+            marker='s',
+            linestyle='none',
+            markersize=15,
+            fillstyle='none'
+        )
+
+        ax.set_ylim(1, 1e6)
+        ax.set_xlim(-1.5, 2.5)
+
+        pdf.savefig(fig)
+        plt.close()
+
+        fig = plt.figure(figsize=(8, 6), tight_layout=True)
+        ax = plt.gca()
+        ax.set_xlabel(r'log$_{10}$($E$/TeV)')
+        ax.set_ylabel(r'$\sigma$ / A')
+
+        ax.errorbar(
+            lge0,
+            rel_err0,
+            marker='o',
+            linestyle='none'
+        )
+        ax.errorbar(
+            lge1,
+            rel_err1,
+            marker='s',
+            linestyle='none',
+            markersize=15,
+            fillstyle='none'
+        )
+
+        ax.set_ylim(0, 0.5)
+        ax.set_xlim(-1.5, 2.5)
+
+        pdf.savefig(fig)
+        plt.close()
+
+        fig = plt.figure(figsize=(8, 6), tight_layout=True)
+        ax = plt.gca()
+        ax.set_xlabel(r'log$_{10}$($E$/TeV)')
+        ax.set_ylabel(r'$\Delta/\sigma$')
+
+        ax.errorbar(
+            lge0,
+            diff,
+            marker='o',
+            linestyle='none'
+        )
+
+        ax.plot(
+            [-1.5, 2.5],
+            [0, 0],
+            color='k',
+            linestyle='--'
+        )
+
+        ax.set_ylim(-0.5, 0.5)
+        ax.set_xlim(-1.5, 2.5)
+
+        pdf.savefig(fig)
+        plt.close()

compressCARE

@@ -93,7 +93,7 @@ merged_care_file=$(merged_care_file "${zenith}" "${atm}" "${wobble}" "${nsb}" "$
 # tar -czvf "${compressed_care_file}" "${merged_dir}"
 # bzip2 -z "${merged_care_file}"
 
-zstd -z "${merged_care_file}"
+zstd -zf "${merged_care_file}"
 
 print_runtime "${SECONDS}"
 

config/GrOptics/GrOpticsV6pilot_generic.txt

@@ -73,3 +73,6 @@ stops with the root prompt. These are cool drawings.
  DRAWTEL 1
 
  TESTTEL 1 PSFgraph
+
+* DEBUGBRANCHES 1
+

countGrOptics

@@ -22,6 +22,16 @@ else
     wobble_list="${wobble}"
 fi
 
+printf "CORSIKAs \n"
+corsika_dir=$(corsika_directory_zip "${zenith}" "${atm}")
+if [ -d "${corsika_dir}" ]; then
+    # shellcheck disable=SC2012
+    n=$(ls "${corsika_dir}" | wc -l)
+    printf "%s" "$n"
+    printf "\n"
+fi        
+
+
 run=$(compute_run "${zenith}" "10")
 
 for wobble in $wobble_list; do

countMerged

@@ -10,9 +10,10 @@ atm="winter"
 zenith="20"
 mode="std"
 wobble=""
+nsb=""
 
 # shellcheck disable=SC2068
-collect_arguments 4 atm zenith mode wobble $@
+collect_arguments 5 atm zenith mode wobble nsb $@
 
 atm=$(validate_atm "${atm}")
 zenith=$(validate_zenith "${zenith}")
@@ -31,7 +32,7 @@ for wobble in ${wobble_list}; do
     merged_care_file=$(merged_care_file "${zenith}" "${atm}" "${wob}" "100" "${mode}")
     merged_care_dir=$(dirname "${merged_care_file}")
 
-    ls "${merged_care_dir}/"*"${wob}"* -lhS
-    n=$(ls "${merged_care_dir}/"*"${wob}"* | wc -l)
+    ls "${merged_care_dir}/"*"${wob}"*"_${nsb}"* -lhS
+    n=$(ls "${merged_care_dir}/"*"${wob}"*"_${nsb}"* | wc -l)
     printf "\n %s files\n" "${n}"
 done

loadFITS.sh

@@ -0,0 +1,8 @@
+#!/bin/zsh
+
+# export VBFSYS=/afs/ifh.de/group/cta/scratch/prado/sw/libVBF
+export FITSSYS=/afs/ifh.de/group/cta/VERITAS/software/FITS/cfitsio
+
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$VBFSYS/lib
+export PATH=$VBFSYS/bin:$VBFSYS/include:$PATH
+