diff --git a/.gitignore b/.gitignore index e43b0f9..74d54ac 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ .DS_Store +tmpdir diff --git a/environment.yml b/environment.yml index 161c0e2..8cf2fae 100644 --- a/environment.yml +++ b/environment.yml @@ -2,6 +2,9 @@ name: eventdisplay-release-tests channels: - conda-forge dependencies: + - python>=3.8 + - numpy + - matplotlib - pandoc - pre-commit - towncrier @@ -11,3 +14,4 @@ dependencies: # activate: conda activate eventdisplay-release-tests # update (conda/mamba): conda env update -f environment.yml --prune # update (micromamba): micromamba update -f environment.yml +# update (micromamba): micromamba update -n eventdisplay-release-tests -f environment.yml diff --git a/release_tests/montecarlo/irf_plotting/irf_plotting.sh b/release_tests/montecarlo/irf_plotting/irf_plotting.sh index 02ddf03..7a41c35 100755 --- a/release_tests/montecarlo/irf_plotting/irf_plotting.sh +++ b/release_tests/montecarlo/irf_plotting/irf_plotting.sh @@ -18,6 +18,7 @@ CUT="NTel3-PointSource-Hard-TMVA-BDT" CUT="NTel2-PointSource-Soft-TMVA-BDT" CUT="NTel2-PointSource-Moderate-TMVA-BDT" CUT="NTel2-PointSource-Moderate" +echo "WARNING: CUT is hardwired to $CUT (for now)" # Comparison plots - version and simtype hardwired COMPAREVERSION="v492" COMPARESIMTYPE="CARE_202404" diff --git a/release_tests/montecarlo/mc_data_comparison/.gitignore b/release_tests/montecarlo/mc_data_comparison/.gitignore deleted file mode 100644 index 15c6bf7..0000000 --- a/release_tests/montecarlo/mc_data_comparison/.gitignore +++ /dev/null @@ -1 +0,0 @@ -tmpdir diff --git a/release_tests/montecarlo/mc_data_comparison/README.md b/release_tests/montecarlo/mc_data_comparison/README.md index a89720b..6a61489 100644 --- a/release_tests/montecarlo/mc_data_comparison/README.md +++ b/release_tests/montecarlo/mc_data_comparison/README.md @@ -1,20 +1,97 @@ -# MC / data comparison +# MC / Data Comparison -Scripts and tools to compare MC distributions with results from the Crab Nebula. +**No apptainer usage!** -Fill and plot distributions with +Scripts and tools to compare Monte Carlo (MC) distributions with observational results from the Crab Nebula. + +## Quick Start ```console -./compareDatawithMC.sh +./compareDatawithMC.sh ``` -Requires as input: +Where `` specifies the observation conditions: + +- `SZE`: Small zenith angles (0-25°, 0.5° wobble) +- `MZE`: Medium zenith angles (40-50°, 0.5° wobble) +- `LZE`: Large zenith angles (50-70°, 0.5° wobble) +- `WOBBLE`: Large wobble offsets (MC at 1° wobble, summer atmosphere only) + +NOTE! The zenith angle bins are different compared to those used in the Crab Nebula analysis. + +## Prerequisites + +### Input Data Requirements + +- MC simulation files for each minor epoch (organized by atmosphere type) +- Crab Nebula mscw results for each minor epoch +- Run lists for Crab observations (generated separately) + +## Preparation Steps + +Before running the comparison, you must generate the Crab run lists: + +```console +# From release_tests/sources/Crab/ +./runlist_generator_from_anasum_log.sh +``` + +This generates run lists for minor epochs, zenith angle ranges, and atmosphere types. +The run lists are saved in `EventDisplay_Release_/Crab/runlists/`. + +## Running the Comparison + +### Step 1: Prepare Run Parameter File + +The `` must contain the following fields (one per line): -- MC files for each minor epoch -- Crab mscw results for each minor epoch +```text +* VERSION # e.g., V6 +* SIMTYPE # e.g., CARE, CARE_RedHV +* ATMOSPHERE # e.g., 61*, 62* +* EPOCH # e.g., V4*, V5*, V6* +* CRAB_NSB # e.g., 0, 1, or NOTSET +``` + +Fields marked with `*` support wildcards for multiple values. + +### Step 2: Execute Comparison + +```console +cd release_tests/montecarlo/mc_data_comparison/ +./compareDatawithMC.sh SZE +``` + +This will: + +1. Parse the run parameter file +2. Locate MC simulation files based on epoch and atmosphere +3. Find corresponding Crab observational data +4. Submit Condor jobs for each combination +5. Generate comparison plots as PDF files + +### Step 3: Monitor Progress + +- Condor jobs are submitted with resource requests: 4000MB RAM, 10GB disk +- Temporary symbolic links are created in: `{output_dir}/tmp/{uuid}/` +- **Note**: Temporary directories must be cleaned up manually after completion + +## Output Structure + +Results are written to: + +```text +EventDisplay_Release_/mc_data_comparison/// +├── _ATM___/ +│ ├── mcdatacomparison.runparameter +│ ├── mcdatacomparison.root +│ ├── mcdatacomparison.log +│ └── *.pdf # Comparison plots +└── tmp/ # Temporary symlinks (manual cleanup required) +``` -**important:** -first run linking of epochs from `../../sources/Crab/`: `./runlist_generator.sh V6` -(or for any other epoch) +## Special Cases -Output and plots are written as PDFs into the `../../../../EventDisplay_Release_/mc_data_comparison/` directories +- **CARE_RedHV simulations**: Handled with modified BDT cut (BDT=0 instead of 1) +- **Summer vs Winter Atmosphere**: Automatically selected based on epoch suffix ('s' or 'w') +- **Wobble Mode**: Forces atmosphere 61 (summer) and 1.0° wobble offset diff --git a/release_tests/montecarlo/mc_data_comparison/compareDatawithMC.sh b/release_tests/montecarlo/mc_data_comparison/compareDatawithMC.sh index 440978b..0a58bca 100755 --- a/release_tests/montecarlo/mc_data_comparison/compareDatawithMC.sh +++ b/release_tests/montecarlo/mc_data_comparison/compareDatawithMC.sh @@ -4,13 +4,13 @@ # # requires: # - MC files for each minor epoch -# - Crab results for each minor epoch (read from Crab runlists; used pre-processed data) +# - Crab results for each minor epoch (read from Crab run lists; used pre-processed data) set -e if [[ $# < 2 ]]; then echo " - ./compareDatawithMC.sh + ./compareDatawithMC.sh --> choose zenith angle / wobble range SZE: small zenith angles (0.5 deg wobble) MZE: medium large zenith angles (0.5 deg wobble) @@ -23,9 +23,9 @@ exit fi ########################### -# read runparameter file +# read run parameter file if [[ ! -e ${1} ]]; then - echo "Error, runparameter file not found: ${1}" + echo "Error, run parameter file not found: ${1}" exit fi # Eventdisplay version @@ -58,10 +58,8 @@ fi # Directory for simulations SIMDIR=${VERITAS_IRFPRODUCTION_DIR}/${VERSION}/${ANALYSISTYPE}/$SIMTYPE/ if [[ ! -e ${SIMDIR} ]]; then - if [[ ! -e ${SIMDIR} ]]; then - echo "Error: simulation directory not found: $SIMDIR" - exit - fi + echo "Error: simulation directory not found: $SIMDIR" + exit fi # Directory for mscw data files DDIR="$VERITAS_PREPROCESSED_DATA_DIR/${ANALYSISTYPE}/mscw/" @@ -123,7 +121,7 @@ do ZEMAX="50." elif [[ $ELE = "LZE" ]] then - ZEMIN="50" + ZEMIN="50." ZEMAX="70." simfile="55deg_${MCWOFF}wob_NOISE${NSB}.mscw.root" # wobble set (everything not 0.5 deg) @@ -216,8 +214,5 @@ do $EVNDISPSCRIPTS/helper_scripts/UTILITY.condorSubmission.sh ${FSCRIPT}.sh 4000M 10G condor_submit ${FSCRIPT}.sh.condor - continue - - done done diff --git a/release_tests/sources/Crab/README.md b/release_tests/sources/Crab/README.md index 34f56af..b1c9330 100644 --- a/release_tests/sources/Crab/README.md +++ b/release_tests/sources/Crab/README.md @@ -43,7 +43,7 @@ for Combine files using pre-processed anasum files and run list generated in step before: ```bash -./anasum_yearly.sh +./anasum_from_runlists.sh ``` ## Plotting diff --git a/release_tests/sources/Crab/plot_all.sh b/release_tests/sources/Crab/plot_all.sh index 5486d36..8564882 100755 --- a/release_tests/sources/Crab/plot_all.sh +++ b/release_tests/sources/Crab/plot_all.sh @@ -19,8 +19,6 @@ ODIR="${2}/$CUT" echo "Output directory: $ODIR" mkdir -p "$ODIR" -PDIR=$(pwd) - # copy anasum log files to release test directory for F in $LFIL; do TF=$(basename $F .combined.root) diff --git a/release_tests/sources/Crab/plot_energy_spectra.C b/release_tests/sources/Crab/plot_energy_spectra.C index 3fc3e6d..882d484 100644 --- a/release_tests/sources/Crab/plot_energy_spectra.C +++ b/release_tests/sources/Crab/plot_energy_spectra.C @@ -14,6 +14,8 @@ #include "TF1.h" +R__LOAD_LIBRARY($EVNDISPSYS/lib/libVAnaSum.so) + #include "../../utilities/parameters.C" #include "../../utilities/printutilities.C" diff --git a/release_tests/sources/Crab/plot_lightcurves.C b/release_tests/sources/Crab/plot_lightcurves.C index 458f068..5b341a7 100644 --- a/release_tests/sources/Crab/plot_lightcurves.C +++ b/release_tests/sources/Crab/plot_lightcurves.C @@ -13,6 +13,8 @@ #include #include +R__LOAD_LIBRARY($EVNDISPSYS/lib/libVAnaSum.so) + #include "../../utilities/parameters.C" #include "../../utilities/printutilities.C" diff --git a/release_tests/sources/Crab/plot_sensitivity.C b/release_tests/sources/Crab/plot_sensitivity.C index b8c4fd8..84de9dd 100644 --- a/release_tests/sources/Crab/plot_sensitivity.C +++ b/release_tests/sources/Crab/plot_sensitivity.C @@ -4,6 +4,7 @@ */ #include +R__LOAD_LIBRARY($EVNDISPSYS/lib/libVAnaSum.so) #include "../../utilities/parameters.C" #include "../../utilities/printutilities.C" diff --git a/release_tests/sources/Crab/plot_sensitivity_from_anasum_log.py b/release_tests/sources/Crab/plot_sensitivity_from_anasum_log.py new file mode 100644 index 0000000..1409d6a --- /dev/null +++ b/release_tests/sources/Crab/plot_sensitivity_from_anasum_log.py @@ -0,0 +1,187 @@ +"""Plot sensitivity diagnostics from anasum log files. + +Reads anasum log files from an input directory and writes a PNG summary figure. + +Usage: python plot_sensitivity_from_anasum_log.py [input_dir] [output_png] +""" + +import argparse +import glob +import os +import re + +import matplotlib.pyplot as plt +import numpy as np + +# Parse input directory argument +parser = argparse.ArgumentParser(description="Plot sensitivity metrics from anasum log files.") +parser.add_argument( + "input_dir", + nargs="?", + default=".", + help="Directory containing anasum log files (default: current directory)", +) +parser.add_argument( + "output_png", + nargs="?", + default="anasum_analysis_plots.png", + help="Output PNG filename (default: anasum_analysis_plots.png)", +) +args = parser.parse_args() + +if not os.path.isdir(args.input_dir): + raise SystemExit(f"Input directory does not exist: {args.input_dir}") + +log_files = sorted( + glob.glob(os.path.join(args.input_dir, "*anasum*V6_20*ATM*SZE_0.5deg.combined.log")) +) + +if not log_files: + raise SystemExit(f"No anasum log files found in: {args.input_dir}") + +epochs = [] +on_off_summary_rates = [] +off_summary_rates = [] +crab_001_times = [] +all_elevations = [] +all_on_rates = [] +all_off_rates = [] + +for filepath in log_files: + print(f"Processing file: {filepath}") + with open(filepath, "r") as f: + content = f.read() + + # Extract epoch from filename (e.g., V6_2023_2023s_ATM62 or V6_2014_2015_ATM61) + basename = os.path.basename(filepath) + epoch_match = re.search(r"(V6_[^_]+_\d+[a-z]*)_ATM\d+", basename) + if not epoch_match: + # Try alternative pattern for ATM61 files + epoch_match = re.search( + r"(V6_[^_]+_\d+[a-z]*)_ATM\d+", basename.replace("_SZE_0.5deg.combined", "") + ) + if not epoch_match: + continue + epoch = epoch_match.group(1) + + # Extract summary on/off rates + rate_match = re.search(r"Rates \(on/Off\): (\d+\.\d+)\s+(\d+\.\d+)", content) + + # Extract 0.01 Crab time in hours from the sensitivity table + lines = content.split("\n") + crab_time = None + in_table = False + for line in lines: + if "[h]" in line: + in_table = True + continue + if in_table and line.strip().startswith("0.01"): + parts = line.split() + if len(parts) >= 3: + try: + crab_time = float(parts[-1]) + except ValueError: + pass + break + + # Extract all elevation and rate values from RUN lines + el_values = [] + on_rates = [] + off_rates = [] + for line in lines: + el_match = re.search(r"RUN\s+\d+.*at\s+(\d+),\s+(-?\d+)\s+deg\s+El\.,\s+Az", line) + rate_match_line = re.search(r"Rates:\s+([\d.]+).*background:\s+([\d.]+)", line) + if el_match: + el_values.append(int(el_match.group(1))) + if rate_match_line: + on_rates.append(float(rate_match_line.group(1))) + off_rates.append(float(rate_match_line.group(2))) + + # Only add data if we have all values and non-empty arrays + if rate_match and crab_time is not None and el_values and on_rates and off_rates: + epochs.append(epoch) + on_off_summary_rates.append(float(rate_match.group(1))) + off_summary_rates.append(float(rate_match.group(2))) + crab_001_times.append(crab_time) + all_elevations.append(np.array(el_values)) + all_on_rates.append(np.array(on_rates)) + all_off_rates.append(np.array(off_rates)) + +# Calculate average elevation per epoch for sensitivity plot +avg_elevations = [np.mean(el) for el in all_elevations] + +# Create plots - 4 rows: 2x2 + 2 extra +fig, axs = plt.subplots(4, 2, figsize=(14, 20)) + +# Plot 1: On Rate vs Epoch +axs[0, 0].plot(epochs, on_off_summary_rates, "o-", color="blue") +axs[0, 0].set_xlabel("Epoch") +axs[0, 0].set_ylabel("On Rate") +axs[0, 0].set_title("On Rate vs Epoch") +axs[0, 0].tick_params(axis="x", rotation=45) + +# Plot 2: Off Rate vs Epoch +axs[0, 1].plot(epochs, off_summary_rates, "s-", color="orange") +axs[0, 1].set_xlabel("Epoch") +axs[0, 1].set_ylabel("Off Rate") +axs[0, 1].set_title("Off Rate vs Epoch") +axs[0, 1].tick_params(axis="x", rotation=45) + +# Plot 3: Violin plot of Elevations per Epoch +axs[1, 0].violinplot(all_elevations, positions=range(len(epochs)), widths=0.8, showmeans=True) +axs[1, 0].set_xticks(range(len(epochs))) +axs[1, 0].set_xticklabels(epochs, rotation=45, ha="right") +axs[1, 0].set_xlabel("Epoch") +axs[1, 0].set_ylabel("Elevation [deg]") +axs[1, 0].set_title("Elevation Distribution per Epoch") + +# Plot 4: Violin plot of On Rates per Epoch +axs[1, 1].violinplot(all_on_rates, positions=range(len(epochs)), widths=0.8, showmeans=True) +axs[1, 1].set_xticks(range(len(epochs))) +axs[1, 1].set_xticklabels(epochs, rotation=45, ha="right") +axs[1, 1].set_xlabel("Epoch") +axs[1, 1].set_ylabel("On Rate") +axs[1, 1].set_title("On Rate Distribution per Epoch") + +# Plot 5: Violin plot of Off Rates per Epoch +axs[2, 0].violinplot(all_off_rates, positions=range(len(epochs)), widths=0.8, showmeans=True) +axs[2, 0].set_xticks(range(len(epochs))) +axs[2, 0].set_xticklabels(epochs, rotation=45, ha="right") +axs[2, 0].set_xlabel("Epoch") +axs[2, 0].set_ylabel("Off Rate") +axs[2, 0].set_title("Off Rate Distribution per Epoch") + +# Plot 6: 0.01 Crab Sensitivity Time vs Epoch +axs[2, 1].plot(epochs, crab_001_times, "s-", color="purple") +axs[2, 1].set_xlabel("Epoch") +axs[2, 1].set_ylabel("0.01 Crab Time [h]") +axs[2, 1].set_title("0.01 Crab Sensitivity Time vs Epoch") +axs[2, 1].tick_params(axis="x", rotation=45) + +# Plot 7: Sensitivity vs Average Elevation +# Color points red if epoch year >= 2020 +colors = [ + ( + "red" + if "2020" in epoch + or "2021" in epoch + or "2022" in epoch + or "2023" in epoch + or "2024" in epoch + or "2025" in epoch + or "2026" in epoch + else "blue" + ) + for epoch in epochs +] +axs[3, 0].scatter(avg_elevations, crab_001_times, c=colors, alpha=0.7) +axs[3, 0].set_xlabel("Average Elevation [deg]") +axs[3, 0].set_ylabel("0.01 Crab Time [h]") +axs[3, 0].set_title("Sensitivity vs Elevation") + +# Hide the empty subplot +axs[3, 1].axis("off") + +plt.tight_layout() +plt.savefig(args.output_png, dpi=300, bbox_inches="tight") +print(f"Processed {len(epochs)} files. Plots saved as {args.output_png}") diff --git a/release_tests/sources/Crab/plot_skymaps.C b/release_tests/sources/Crab/plot_skymaps.C index 98d905a..5c8fda3 100644 --- a/release_tests/sources/Crab/plot_skymaps.C +++ b/release_tests/sources/Crab/plot_skymaps.C @@ -11,6 +11,8 @@ #include #include +R__LOAD_LIBRARY($EVNDISPSYS/lib/libVAnaSum.so) + #include "../../utilities/parameters.C" #include "../../utilities/printutilities.C" diff --git a/release_tests/sources/Crab/runlist_generator_from_anasum_log.sh b/release_tests/sources/Crab/runlist_generator_from_anasum_log.sh index 913b3b1..7783734 100755 --- a/release_tests/sources/Crab/runlist_generator_from_anasum_log.sh +++ b/release_tests/sources/Crab/runlist_generator_from_anasum_log.sh @@ -99,6 +99,11 @@ do echo "Run $R - log file not found: ${ANASUMLOG}" exit fi + # Skip runs without 4-telescope cuts + if ! grep -q "VGammaHadronCuts::printCutSummary() (ntel=4" "${ANASUMLOG}"; then + echo "RUN $R not a 4-telescope run..skipping" + continue + fi echo "DATADIR ${ANASUMLOG}" # read and extract run info from files INSTRUMENT_EPOCH=$(grep "Instrument epoch selected" "${ANASUMLOG}" | head -n 1 | awk '{print $NF}') @@ -117,7 +122,7 @@ do OBSL="stdHV" fi ELEV=$(grep "mean elevation" "${ANASUMLOG}" | head -n 1 | awk '{print $3}') - EL=$(echo $ELEV | awk -v e=$ELEV '{if (e > 50 ) {print "SZE"} else if (e > 40 ) {print "MZE"} else if (e > 30 ) {print "LZE"} else {print "BZE"}}') + EL=$(echo $ELEV | awk -v e=$ELEV '{if (e > 60 ) {print "SZE"} else if (e > 45 ) {print "MZE"} else if (e > 35 ) {print "LZE"} else {print "BZE"}}') WOBBLESTRING=$(grep "Wobble offsets (currE)" "${ANASUMLOG}") n_offset=$(echo "$WOBBLESTRING" | head -n 1 | awk '{print $5}') w_offset=$(echo "$WOBBLESTRING" | head -n 1 | awk '{print $7}') diff --git a/release_tests/sources/Crab/runlist_releaseTestingV6.dat b/release_tests/sources/Crab/runlist_releaseTestingV6.dat index da6d0e3..e5d6a73 100644 --- a/release_tests/sources/Crab/runlist_releaseTestingV6.dat +++ b/release_tests/sources/Crab/runlist_releaseTestingV6.dat @@ -814,7 +814,6 @@ 104174 104199 104201 -104202 104251 104252 104253 @@ -992,7 +991,6 @@ 113299 113303 113352 -113375 113515 113637 113689 diff --git a/release_tests/utilities/parameters.C b/release_tests/utilities/parameters.C index a6880a8..cf49eeb 100644 --- a/release_tests/utilities/parameters.C +++ b/release_tests/utilities/parameters.C @@ -17,8 +17,9 @@ * Helper for ROOT macros: load libVAnaSum from environment. * Search order: * 1) $VERITAS_VANASUM_LIBRARY - * 2) $EVNDISP/lib/libVAnaSum.so - * 3) ROOT library path via "libVAnaSum.so" + * 2) $EVNDISPSYS/lib/libVAnaSum.so + * 3) $EVNDISP/lib/libVAnaSum.so + * 4) ROOT library path via "libVAnaSum.so" */ bool loadVAnaSumLibrary() { @@ -35,6 +36,14 @@ bool loadVAnaSumLibrary() iLibPath = iEnvLib; } if( iLibPath.size() == 0 ) + { + const char* iEvndispSys = gSystem->Getenv( "EVNDISPSYS" ); + if( iEvndispSys ) + { + iLibPath = string( iEvndispSys ) + "/lib/libVAnaSum.so"; + } + } + if( iLibPath.size() == 0 ) { const char* iEvndisp = gSystem->Getenv( "EVNDISP" ); if( iEvndisp ) @@ -45,21 +54,26 @@ bool loadVAnaSumLibrary() if( iLibPath.size() > 0 && !gSystem->AccessPathName( iLibPath.c_str() ) ) { - if( gSystem->Load( iLibPath.c_str() ) >= 0 ) + // Use gROOT->ProcessLine to properly load the library and its dictionaries + string loadCmd = "R__LOAD_LIBRARY(" + iLibPath + ");"; + int result = gROOT->ProcessLine( loadCmd.c_str() ); + if( result >= 0 ) { iLibraryLoaded = true; return true; } } - if( gSystem->Load( "libVAnaSum.so" ) >= 0 ) + // Try with just the library name + int result = gROOT->ProcessLine( "R__LOAD_LIBRARY(libVAnaSum.so);" ); + if( result >= 0 ) { iLibraryLoaded = true; return true; } cout << "Error: unable to load libVAnaSum.so. " - << "Set VERITAS_VANASUM_LIBRARY or EVNDISP." << endl; + << "Set VERITAS_VANASUM_LIBRARY, EVNDISPSYS, or EVNDISP." << endl; return false; }