diff --git a/.travis.yml b/.travis.yml index 8a24fa2e..23a63d5e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,8 +11,11 @@ env: - ROOT=5-34-18 - ROOT=6-00-00 global: - - CLANG_VERSION=3.4 + - CLANG_VERSION=3.7 + - GCC_VERSION=4.9 + - BOOST_VERSION=1.55 +before_install: source ci/before_install.sh install: source ci/install.sh script: bash ci/test.sh cache: apt diff --git a/ci/before_install.sh b/ci/before_install.sh new file mode 100644 index 00000000..35f44810 --- /dev/null +++ b/ci/before_install.sh @@ -0,0 +1,41 @@ +#!/bin/bash +# This script is meant to be called by the "before_install" step defined in +# .travis.yml. See http://docs.travis-ci.com/ for more details. +# The behaviour of the script is controlled by environment variabled defined +# in the .travis.yml in the top level folder of the project. + +#set -e + +# Check if we are running Python 2 or 3. This is needed for the apt-get package names +if [[ $TRAVIS_PYTHON_VERSION == '3.2' ]]; then + export PYTHON_SUFFIX="3"; +fi + +sudo add-apt-repository --yes ppa:kalakris/cmake +# add repositories for gcc ${GCC_VERSION} and clang $CLANG_VERSION (set in .travis.yml) +sudo add-apt-repository --yes ppa:ubuntu-toolchain-r/test +sudo add-apt-repository --yes 'deb http://llvm.org/apt/precise/ llvm-toolchain-precise main' +sudo add-apt-repository --yes 'deb http://ppa.launchpad.net/boost-latest/ppa/ubuntu precise main' +wget -O - http://llvm.org/apt/llvm-snapshot.gpg.key | sudo apt-key add - +# Needed because sometimes travis' repositories get out of date +sudo apt-get update -q +# Install the dependencies we need +time sudo apt-get -q install cmake clang-${CLANG_VERSION} libclang-${CLANG_VERSION}-dev gcc-${GCC_VERSION} g++-${GCC_VERSION} boost${BOOST_VERSION} libboost${BOOST_VERSION}-dev libboost-test${BOOST_VERSION}-dev \ +python${PYTHON_SUFFIX}-numpy python${PYTHON_SUFFIX}-sphinx python${PYTHON_SUFFIX}-nose python${PYTHON_SUFFIX}-pip cython${PYTHON_SUFFIX} +if [[ $TRAVIS_PYTHON_VERSION == '2.7' ]]; then + time sudo apt-get install -qq python-matplotlib python-tables; +fi + +pip install nose --upgrade +nosetests -V + +# setup newer compilers ( we need gcc >= 4.7 for c++11 +sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-${GCC_VERSION} 50; +sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-${GCC_VERSION} 50; +sudo update-alternatives --set gcc /usr/bin/gcc-${GCC_VERSION}; +sudo update-alternatives --set g++ /usr/bin/g++-${GCC_VERSION}; + +sudo update-alternatives --install /usr/bin/clang clang /usr/bin/clang-${CLANG_VERSION} 50; +sudo update-alternatives --install /usr/bin/clang++ clang++ /usr/bin/clang++-${CLANG_VERSION} 50; +sudo update-alternatives --set clang /usr/bin/clang-${CLANG_VERSION}; +sudo update-alternatives --set clang++ /usr/bin/clang++-${CLANG_VERSION}; diff --git a/ci/install.sh b/ci/install.sh index c84abf7e..316eb23d 100755 --- a/ci/install.sh +++ b/ci/install.sh @@ -4,27 +4,7 @@ # The behavior of the script is controlled by environment variabled defined # in the .travis.yml in the top level folder of the project. -#set -e - -# Check if we are running Python 2 or 3. This is needed for the apt-get package names -if [[ $TRAVIS_PYTHON_VERSION == '3.2' ]]; then - export PYTHON_SUFFIX="3"; -fi - -# add repositories for gcc 4.8 and clang $CLANG_VERSION (set in .travis.yml) -sudo add-apt-repository --yes ppa:ubuntu-toolchain-r/test -sudo add-apt-repository --yes 'deb http://llvm.org/apt/precise/ llvm-toolchain-precise main' -wget -O - http://llvm.org/apt/llvm-snapshot.gpg.key | sudo apt-key add - -# Needed because sometimes travis' repositories get out of date -time sudo apt-get update -qq - -# Install the dependencies we need -time sudo apt-get -qq install clang-${CLANG_VERSION} libclang-${CLANG_VERSION}-dev gcc-4.8 g++-4.8 -time sudo apt-get install -qq python${PYTHON_SUFFIX}-numpy python${PYTHON_SUFFIX}-sphinx python${PYTHON_SUFFIX}-nose python${PYTHON_SUFFIX}-pip cython${PYTHON_SUFFIX} -# matplotlib and PyTables are not available for Python 3 as packages from the main repo yet. -if [[ $TRAVIS_PYTHON_VERSION == '2.7' ]]; then - time sudo apt-get install -qq python-matplotlib python-tables; -fi +set -e # Install a ROOT binary that we custom-built in a 64-bit Ubuntu VM # for the correct python / ROOT version @@ -38,38 +18,16 @@ source root/bin/thisroot.sh root -l -q python -c "import ROOT; ROOT.TBrowser()" -# setup newer compilers for ROOT 6 -if [[ $ROOT == '6-00-00' ]]; then - sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-4.8 50; - sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-4.8 50; - sudo update-alternatives --set gcc /usr/bin/gcc-4.8; sudo update-alternatives --set g++ /usr/bin/g++-4.8; -fi - -# setup vpython with all packages -# being standalone.sh #DailyPythonTools location export base=`pwd` +# setup vpython with all packages # package list from FinalStateAnalysis # (https://github.com/uwcms/FinalStateAnalysis/blob/master/recipe/install_python.sh) -echo "Installing yolk" -time sudo pip install -U yolk -echo "Installing ipython" -time sudo pip install -U ipython -echo "Installing termcolor" -time sudo pip install -U termcolor echo "Installing uncertainties <-- awesome error propagation" -time sudo pip install -U uncertainties -echo "Install progressbar" -time sudo pip install -U progressbar -echo "Installing argparse" -time sudo pip install -U argparse -echo "Installing pudb <-- interactive debugging" -time sudo pip install -U pudb -echo "Installing dateutil" -time sudo pip install python-dateutil -echo "Installing PrettyTable" -time sudo pip install PrettyTable +time pip install -U uncertainties +echo "Installing tabulate (latex table printing, etc)" +pip install tabulate echo "Installing rootpy" time pip install --user -e $base/external/rootpy @@ -82,7 +40,8 @@ if [ ! -d "$base/external/lib" ]; then mkdir $base/external/lib echo "Building RooUnfold" cd $base/external/RooUnfold/ - make -j4 + cmake CMakeLists.txt + make #remove tmp folder rm -fr $base/external/RooUnfold/tmp mv $base/external/RooUnfold/libRooUnfold.so $base/external/lib/. @@ -92,7 +51,7 @@ if [ ! -d "$base/external/lib" ]; then echo "Building TopAnalysis" cd $base/external/TopAnalysis/ - make -j4 + make # remove tmp folder rm -fr $base/external/TopAnalysis/tmp mv $base/external/TopAnalysis/libTopSVDUnfold.so $base/external/lib/. @@ -102,6 +61,6 @@ fi cd $base export PATH=$PATH:$base/bin -# end standalone.sh + # add base path from setup_standalone to PYTHONPATH export PYTHONPATH=$PYTHONPATH:$base diff --git a/config/unfolding/HT_7TeV_electron_channel.json b/config/unfolding/HT_7TeV_electron_channel.json new file mode 100644 index 00000000..123afd45 --- /dev/null +++ b/config/unfolding/HT_7TeV_electron_channel.json @@ -0,0 +1,24 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_electron_channel/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_electron_channel/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_electron_channel/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/HT/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "electron", + "variable": "HT" +} \ No newline at end of file diff --git a/config/unfolding/HT_7TeV_muon_channel.json b/config/unfolding/HT_7TeV_muon_channel.json new file mode 100644 index 00000000..52c18a5a --- /dev/null +++ b/config/unfolding/HT_7TeV_muon_channel.json @@ -0,0 +1,24 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_muon_channel/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_muon_channel/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_muon_channel/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/HT/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "muon", + "variable": "HT" +} \ No newline at end of file diff --git a/config/unfolding/HT_8TeV_electron_channel.json b/config/unfolding/HT_8TeV_electron_channel.json new file mode 100644 index 00000000..e34b426e --- /dev/null +++ b/config/unfolding/HT_8TeV_electron_channel.json @@ -0,0 +1,24 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_electron_channel/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_electron_channel/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_electron_channel/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/HT/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "electron", + "variable": "HT" +} \ No newline at end of file diff --git a/config/unfolding/HT_8TeV_muon_channel.json b/config/unfolding/HT_8TeV_muon_channel.json new file mode 100644 index 00000000..9dec4bf0 --- /dev/null +++ b/config/unfolding/HT_8TeV_muon_channel.json @@ -0,0 +1,24 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_muon_channel/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_muon_channel/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_HT_analyser_muon_channel/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/HT/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "muon", + "variable": "HT" +} \ No newline at end of file diff --git a/config/unfolding/MET_7TeV_electron_channel.json b/config/unfolding/MET_7TeV_electron_channel.json new file mode 100644 index 00000000..ba8c17fc --- /dev/null +++ b/config/unfolding/MET_7TeV_electron_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/MET/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "electron", + "variable": "MET", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/MET_7TeV_muon_channel.json b/config/unfolding/MET_7TeV_muon_channel.json new file mode 100644 index 00000000..4801e82d --- /dev/null +++ b/config/unfolding/MET_7TeV_muon_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_muon_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_muon_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_muon_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/MET/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "muon", + "variable": "MET", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/MET_8TeV_electron_channel.json b/config/unfolding/MET_8TeV_electron_channel.json new file mode 100644 index 00000000..b34c60fd --- /dev/null +++ b/config/unfolding/MET_8TeV_electron_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/MET/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "electron", + "variable": "MET", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/MET_8TeV_muon_channel.json b/config/unfolding/MET_8TeV_muon_channel.json new file mode 100644 index 00000000..81d81ffc --- /dev/null +++ b/config/unfolding/MET_8TeV_muon_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_muon_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_muon_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MET_analyser_muon_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/MET/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "muon", + "variable": "MET", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/MT_7TeV_electron_channel.json b/config/unfolding/MT_7TeV_electron_channel.json new file mode 100644 index 00000000..2fd884d3 --- /dev/null +++ b/config/unfolding/MT_7TeV_electron_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/MT/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "electron", + "variable": "MT", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/MT_7TeV_muon_channel.json b/config/unfolding/MT_7TeV_muon_channel.json new file mode 100644 index 00000000..b07d5dd9 --- /dev/null +++ b/config/unfolding/MT_7TeV_muon_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_muon_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_muon_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_muon_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/MT/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "muon", + "variable": "MT", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/MT_8TeV_electron_channel.json b/config/unfolding/MT_8TeV_electron_channel.json new file mode 100644 index 00000000..e55575ee --- /dev/null +++ b/config/unfolding/MT_8TeV_electron_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/MT/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "electron", + "variable": "MT", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/MT_8TeV_muon_channel.json b/config/unfolding/MT_8TeV_muon_channel.json new file mode 100644 index 00000000..6a240d43 --- /dev/null +++ b/config/unfolding/MT_8TeV_muon_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_muon_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_muon_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_MT_analyser_muon_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/MT/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "muon", + "variable": "MT", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/ST_7TeV_electron_channel.json b/config/unfolding/ST_7TeV_electron_channel.json new file mode 100644 index 00000000..0ec75811 --- /dev/null +++ b/config/unfolding/ST_7TeV_electron_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/ST/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "electron", + "variable": "ST", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/ST_7TeV_muon_channel.json b/config/unfolding/ST_7TeV_muon_channel.json new file mode 100644 index 00000000..8f11cd74 --- /dev/null +++ b/config/unfolding/ST_7TeV_muon_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_muon_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_muon_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_muon_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/ST/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "muon", + "variable": "ST", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/ST_8TeV_electron_channel.json b/config/unfolding/ST_8TeV_electron_channel.json new file mode 100644 index 00000000..fa05fd27 --- /dev/null +++ b/config/unfolding/ST_8TeV_electron_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/ST/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "electron", + "variable": "ST", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/ST_8TeV_muon_channel.json b/config/unfolding/ST_8TeV_muon_channel.json new file mode 100644 index 00000000..780399f4 --- /dev/null +++ b/config/unfolding/ST_8TeV_muon_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_muon_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_muon_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_ST_analyser_muon_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/ST/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "muon", + "variable": "ST", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/WPT_7TeV_electron_channel.json b/config/unfolding/WPT_7TeV_electron_channel.json new file mode 100644 index 00000000..1487004e --- /dev/null +++ b/config/unfolding/WPT_7TeV_electron_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/WPT/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "electron", + "variable": "WPT", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/WPT_7TeV_muon_channel.json b/config/unfolding/WPT_7TeV_muon_channel.json new file mode 100644 index 00000000..aaf8d6e5 --- /dev/null +++ b/config/unfolding/WPT_7TeV_muon_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_muon_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_muon_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_muon_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/7TeV/WPT/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/7TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 7, + "channel": "muon", + "variable": "WPT", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/WPT_8TeV_electron_channel.json b/config/unfolding/WPT_8TeV_electron_channel.json new file mode 100644 index 00000000..f54b8897 --- /dev/null +++ b/config/unfolding/WPT_8TeV_electron_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/WPT/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "electron", + "variable": "WPT", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/WPT_8TeV_muon_channel.json b/config/unfolding/WPT_8TeV_muon_channel.json new file mode 100644 index 00000000..77ba1df3 --- /dev/null +++ b/config/unfolding/WPT_8TeV_muon_channel.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_muon_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_muon_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root", + "histogram": "unfolding_WPT_analyser_muon_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/WPT/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "muon", + "variable": "WPT", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/systematics/MET_8TeV_electron_channel_scale_down.json b/config/unfolding/systematics/MET_8TeV_electron_channel_scale_down.json new file mode 100644 index 00000000..543cc4ba --- /dev/null +++ b/config/unfolding/systematics/MET_8TeV_electron_channel_scale_down.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_scaledown_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_scaledown_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_scaledown_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/MET/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "electron", + "variable": "MET", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/config/unfolding/systematics/MET_8TeV_electron_channel_scale_up.json b/config/unfolding/systematics/MET_8TeV_electron_channel_scale_up.json new file mode 100644 index 00000000..543cc4ba --- /dev/null +++ b/config/unfolding/systematics/MET_8TeV_electron_channel_scale_up.json @@ -0,0 +1,26 @@ +{ + "config_for" : "src/unfolding_tests/get_best_regularisation.py", + "truth" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_scaledown_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/truth" + }, + "gen_vs_reco" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_scaledown_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/response_without_fakes" + }, + "measured" : { + "file": "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_scaledown_asymmetric.root", + "histogram": "unfolding_MET_analyser_electron_channel_patType1CorrectedPFMet/measured" + }, + "data" : { + "file": "data/absolute_eta_M3_angle_bl/8TeV/MET/fit_results/central/fit_results_electron_patType1CorrectedPFMet.txt", + "histogram": "TTJet" + }, + "output_folder": "plots/8TeV/unfolding_tests", + "output_format": ["png", "pdf"], + "centre-of-mass energy" : 8, + "channel": "electron", + "variable": "MET", + "n_toy": 1000, + "n_tau_scan_points": 100 +} \ No newline at end of file diff --git a/examples/forest/README.md b/examples/forest/README.md new file mode 100644 index 00000000..b561ec94 --- /dev/null +++ b/examples/forest/README.md @@ -0,0 +1,2 @@ +# examples/forest +Everything concerning Trees and NTuples \ No newline at end of file diff --git a/examples/forest/random_tree.py b/examples/forest/random_tree.py new file mode 100644 index 00000000..f7549eef --- /dev/null +++ b/examples/forest/random_tree.py @@ -0,0 +1,23 @@ +from rootpy.tree import Tree +from rootpy.io import root_open +from random import gauss + +f = root_open("test.root", "recreate") + +tree = Tree("test") +tree.create_branches( + {'x': 'F', + 'y': 'F', + 'z': 'F', + 'i': 'I'}) + +for i in xrange(10000): + tree.x = gauss(.5, 1.) + tree.y = gauss(.3, 2.) + tree.z = gauss(13., 42.) + tree.i = i + tree.fill() +tree.write() + +f.close() + diff --git a/examples/read_ntuples.py b/examples/forest/read_ntuples.py similarity index 100% rename from examples/read_ntuples.py rename to examples/forest/read_ntuples.py diff --git a/experimental/README.md b/experimental/README.md new file mode 100644 index 00000000..ba3f427c --- /dev/null +++ b/experimental/README.md @@ -0,0 +1,4 @@ +# experimental +Home to experimental (for incubation) scripts, modules and examples that don't really fit anywhere else. +Please use scripts with caution, this is potentially a mine field. +If a script/module/example is stable enough (reached maturity) it should be move to an appropriate location. diff --git a/experimental/TMatrix.py b/experimental/TMatrix.py new file mode 100644 index 00000000..24b5380a --- /dev/null +++ b/experimental/TMatrix.py @@ -0,0 +1,26 @@ +from ROOT import TMatrixDSparse +from array import array +n_x, n_y = 4,4 +data_array = array('d', [i+1 for i in range (n_x*n_y)]) +row = array('i', [i for i in range (n_x)]) +col = array('i', [i for i in range (n_y)]) +print col, row +m_square = TMatrixDSparse(n_x,n_y) +m_square.SetSparseIndex(n_x) +m_square.SetRowIndexArray(row) +m_square.SetColIndexArray(col) +# m_square.SetMatrixArray(5, row, col, data_array) + +row_index = m_square.GetRowIndexArray() +col_index = m_square.GetColIndexArray() +pdata = m_square.GetMatrixArray() +print 'n rows=', m_square.GetNrows() + +for i in range(m_square.GetNrows()): + sIndex = row_index[i] + eIndex = row_index[i+1] + print sIndex, eIndex + for index in range(sIndex, eIndex): + icol = col_index[index] + data = pdata[index] + print 'data(%d,%d) = %.4f' % (i, icol, data) diff --git a/experimental/W_studies/out.txt b/experimental/W_studies/out.txt new file mode 100644 index 00000000..1ac88117 --- /dev/null +++ b/experimental/W_studies/out.txt @@ -0,0 +1,11 @@ +nohup: ignoring input +[?1034h------ EPlusJets ------ +--- Inclusive +--- PU Binned +--- Leading jet eta +--- Sub leading jet eta +------ MuPlusJets ------ +--- Inclusive +--- PU Binned +--- Leading jet eta +--- Sub leading jet eta diff --git a/experimental/condor/01_fit.description b/experimental/condor/01_fit.description new file mode 100644 index 00000000..6857ba1e --- /dev/null +++ b/experimental/condor/01_fit.description @@ -0,0 +1,18 @@ +Executable = experimental/condor/run_01.sh +Universe = vanilla +Output = 01.job.$(cluster).$(process).out +Error = 01.job.$(cluster).$(process).err +Log = 01.job.$(cluster).$(process).log +arguments = $(process) + +transfer_input_files = dps.tar +should_transfer_files = YES +when_to_transfer_output = ON_EXIT_OR_EVICT + +request_cpus = 1 +request_memory=500 + +# use the ENV that is provided +getenv = true + +queue 4 diff --git a/experimental/condor/run01_forAllOptions.py b/experimental/condor/run01_forAllOptions.py new file mode 100644 index 00000000..c0af8d4a --- /dev/null +++ b/experimental/condor/run01_forAllOptions.py @@ -0,0 +1,23 @@ +import os +from optparse import OptionParser + +fit_var="M3,angle_bl" + +vars = [ + 'MET', + 'HT', + 'ST', + 'WPT' +] + +jobOptions = [ '-v %s --fit-variables \'%s\' ' % ( var, fit_var ) for var in vars ] + +parser = OptionParser("Merge histogram files on DICE") +parser.add_option("-n", dest="jobNumber", default=-1, type='int', + help="Specify which job number to run") + +(options, _) = parser.parse_args() + +jobOption = jobOptions[options.jobNumber] +print 'Running with options : ',jobOption +os.system('python src/cross_section_measurement/01_get_fit_results.py %s' % jobOption ) diff --git a/experimental/condor/run_01.sh b/experimental/condor/run_01.sh new file mode 100644 index 00000000..24ee5025 --- /dev/null +++ b/experimental/condor/run_01.sh @@ -0,0 +1,22 @@ +#!/bin/bash +scramv1 project CMSSW CMSSW_7_4_0_pre7 +cd CMSSW_7_4_0_pre7/src +eval `scramv1 runtime -sh` +tar -xf ../../dps.tar +cd DailyPythonScripts/ +cd external/ +rm -r vpython/ +cd ../ +git submodule init && git submodule update +./setup.sh +eval `scramv1 runtime -sh` +. environment.sh +python experimental/condor/run01_forAllOptions.py -n $1 +echo "Done" +ls +echo "Tarring" +tar -cf output_$1.tar data +ls +echo "Moving" +mv output_$1.tar ../../../ +ls ../../../ \ No newline at end of file diff --git a/experimental/hist_to_numpy.py b/experimental/hist_to_numpy.py new file mode 100644 index 00000000..131ebe38 --- /dev/null +++ b/experimental/hist_to_numpy.py @@ -0,0 +1,31 @@ +''' +Created on 28 Jul 2014 + +@author: phxlk +''' +import numpy as np +from rootpy.plotting import Hist + +def hist_to_numpy(hist): + edges = list(hist.xedges()) + np_hist, be = np.histogram(list(hist.y()), bins = hist.nbins(), range = (edges[0], edges[-1])) + return np_hist, be + + +n_bins = 100 +min_x = 0 +max_x = 200 +N_bkg1_ctl = 30000 +N_signal_ctl = 2000 +N_bkg1_obs = 30000 +N_signal_obs = 2000 +N_data = N_bkg1_obs + N_signal_obs +mu1, mu2, sigma1, sigma2 = 100, 140, 15, 5 +mu3, mu4, sigma3, sigma4 = 80, 170, 14, 10 + +x1 = mu1 + sigma1 * np.random.randn( N_bkg1_ctl ) +h1 = Hist( n_bins, 0, 200, title = 'data' ) +map( h1.Fill, x1 ) + + + diff --git a/experimental/minuit.py b/experimental/minuit.py new file mode 100644 index 00000000..eea0ba27 --- /dev/null +++ b/experimental/minuit.py @@ -0,0 +1,60 @@ +from iminuit import describe +from tools.Fitting import IMinuit, FitData, FitDataCollection +from rootpy.plotting import Hist +import numpy as np + +N_bkg1 = 9000 +N_signal = 1000 +N_bkg1_obs = 10000 +N_signal_obs = 2000 +N_data = N_bkg1_obs + N_signal_obs +mu1, mu2, sigma1, sigma2 = 100, 140, 15, 5 +x1 = mu1 + sigma1 * np.random.randn( N_bkg1 ) +x2 = mu2 + sigma2 * np.random.randn( N_signal ) +x1_obs = mu1 + sigma1 * np.random.randn( N_bkg1_obs ) +x2_obs = mu2 + sigma2 * np.random.randn( N_signal_obs ) + +x3 = mu2 + sigma1 * np.random.randn( N_bkg1 ) +x4 = mu1 + sigma2 * np.random.randn( N_signal ) +x3_obs = mu2 + sigma1 * np.random.randn( N_bkg1_obs ) +x4_obs = mu1 + sigma2 * np.random.randn( N_signal_obs ) + +data_scale = 1.2 +N_data = N_data * data_scale + +h_bkg1_1 = Hist( 100, 40, 200, title = 'Background' ) +h_signal_1 = h_bkg1_1.Clone( title = 'Signal' ) +h_data_1 = h_bkg1_1.Clone( title = 'Data' ) +h_bkg1_2 = h_bkg1_1.Clone( title = 'Background' ) +h_signal_2 = h_bkg1_1.Clone( title = 'Signal' ) +h_data_2 = h_bkg1_1.Clone( title = 'Data' ) + +# fill the histograms with our distributions +map( h_bkg1_1.Fill, x1 ) +map( h_signal_1.Fill, x2 ) +map( h_data_1.Fill, x1_obs ) +map( h_data_1.Fill, x2_obs ) + +map( h_bkg1_2.Fill, x3 ) +map( h_signal_2.Fill, x4 ) +map( h_data_2.Fill, x3_obs ) +map( h_data_2.Fill, x4_obs ) + +h_data_1.Scale( data_scale ) +h_data_2.Scale( data_scale ) + +histograms_1 = {'signal': h_signal_1, + 'bkg1': h_bkg1_1} + +histograms_2 = {'signal': h_signal_2, + 'bkg1': h_bkg1_2} + +fit_data_1 = FitData( h_data_1, histograms_1, fit_boundaries = ( 40, 200 ) ) +fit_data_2 = FitData( h_data_2, histograms_2, fit_boundaries = ( 40, 200 ) ) + +single_fit_collection = FitDataCollection() +single_fit_collection.add( fit_data_1 ) + +m = IMinuit(single_fit_collection) +print describe(m.likelihood_3_samples, verbose=True) +print describe(m.likelihood_4_samples, verbose=True) \ No newline at end of file diff --git a/experimental/probfit/basics.py b/experimental/probfit/basics.py new file mode 100644 index 00000000..6d48052a --- /dev/null +++ b/experimental/probfit/basics.py @@ -0,0 +1,97 @@ +import numpy as np +from numpy import array +from rootpy.plotting import Hist +from iminuit import Minuit +from probfit import BinnedLH, Extended, AddPdf, draw_blh +from probfit.pdf import HistogramPdf +from probfit.plotting import draw_pdf +from root_numpy import root2array + +n_bins = 100 +min_x = 0 +max_x = 200 +N_bkg1_ctl = 30000 +N_signal_ctl = 2000 +N_bkg1_obs = 30000 +N_signal_obs = 2000 +N_data = N_bkg1_obs + N_signal_obs +mu1, mu2, sigma1, sigma2 = 100, 140, 15, 5 +mu3, mu4, sigma3, sigma4 = 80, 170, 14, 10 + +def main(): + data, data_ctl = get_data() + signal, bkg, signal_ctl, bkg_ctl = get_templates() + edges = list( data.xedges() ) +# print list(data.y()) +# data, be = hist_to_numpy(data) +# print data +# data_ctl, _ = hist_to_numpy(data_ctl) +# signal, _ = hist_to_numpy(signal) +# signal_ctl, _ = hist_to_numpy(signal_ctl) +# bkg, _ = hist_to_numpy(bkg) +# bkg_ctl, _ = hist_to_numpy(bkg_ctl) + sig_pdf = hist_to_pdf( signal ) + bkg_pdf = hist_to_pdf( bkg ) + epsig = Extended( sig_pdf, extname = 'N1' ) + epbkg = Extended( bkg_pdf, extname = 'N2' ) + pdf = AddPdf( epbkg, epsig ) + fitdata = array( list( data.y() ) ) + blh = BinnedLH( pdf, fitdata, bins = n_bins, bound = ( min_x, max_x ), extended = True ) + m = Minuit( blh, N1 = 2000, N2 = 30000, error_N1 = 44, error_N2 = 200 ) + m.migrad() + blh.draw( m, parts = True ) + + +def get_data(): + # start with data that has the same statistics + + # randomise it + x1_ctl = mu1 + sigma1 * np.random.randn( N_bkg1_ctl ) + x2_ctl = mu2 + sigma2 * np.random.randn( N_signal_ctl ) + x1_obs = mu3 + sigma3 * np.random.randn( N_bkg1_obs ) + x2_obs = mu4 + sigma4 * np.random.randn( N_signal_obs ) + + h1 = Hist( n_bins, 0, 200, title = 'data' ) + h2 = Hist( n_bins, 0, 200, title = 'data_ctl' ) + + # fill the histograms with our distributions + map( h1.Fill, x1_obs ) + map( h1.Fill, x2_obs ) + map( h2.Fill, x1_ctl ) + map( h2.Fill, x2_ctl ) + return h1, h2 + +def get_templates(): + x1_ctl = mu1 + sigma1 * np.random.randn( N_bkg1_ctl ) + x2_ctl = mu2 + sigma2 * np.random.randn( N_signal_ctl ) + x1_obs = mu3 + sigma3 * np.random.randn( N_bkg1_obs * 2 ) + x2_obs = mu4 + sigma4 * np.random.randn( N_signal_obs * 2 ) + + signal_1 = Hist( n_bins, 0, 200, title = 'signal_1' ) + bkg_1 = Hist( n_bins, 0, 200, title = 'bkg_1' ) + + signal_2 = Hist( n_bins, 0, 200, title = 'signal_1' ) + bkg_2 = Hist( n_bins, 0, 200, title = 'bkg_1' ) + + # fill the histograms with our distributions + map( bkg_1.Fill, x1_obs ) + map( signal_1.Fill, x2_obs ) + map( bkg_2.Fill, x1_ctl ) + map( signal_2.Fill, x2_ctl ) + + + return signal_1, bkg_1, signal_2, bkg_2 + +def hist_to_numpy( hist ): + edges = list( hist.xedges() ) + np_hist, be = np.histogram( list( hist.y() ), bins = hist.nbins(), range = ( edges[0], edges[-1] ) ) + return np_hist, be + +def hist_to_pdf( hist ): + edges = array( list( hist.xedges() ) ) + y_values = array( list( hist.y() ) ) + return HistogramPdf( y_values, edges ) + +if __name__ == "__main__": + main() + diff --git a/experimental/probfit/test.py b/experimental/probfit/test.py new file mode 100644 index 00000000..b4ab54a5 --- /dev/null +++ b/experimental/probfit/test.py @@ -0,0 +1,48 @@ +from iminuit import Minuit +from probfit import BinnedLH, gaussian, Extended +from matplotlib import pyplot as plt +from numpy.random import randn +from rootpy.plotting import Hist +from numpy import array +import numpy as np + +data = randn(1000)*2 + 2 +# print data +h1 = Hist( 40, -6, 10, title = 'data' ) +map( h1.Fill, data ) +# data = array( list( h1.y() ) ) +# print data +#Unextended + +# blh = BinnedLH(gaussian, data, bins = 40, bound = ( -6, 10 )) +# # #if you wonder what it loos like call desceibe(blh) +# m = Minuit(blh, mean=0., sigma=0.5) +# +# plt.figure(figsize=(8, 6)) +# plt.subplot(221) +# blh.draw(m) +# plt.title('Unextended Before') +# +# m.migrad() # fit +# +# plt.subplot(222) +# blh.draw(m) +# plt.title('Unextended After') +# +# #Extended +# +# ext_gauss = Extended(gaussian) +# +# blh = BinnedLH(ext_gauss, data, extended=True, bins = 40, bound = ( -6, 10 )) +# m = Minuit(blh, mean=0., sigma=0.5, N=900.) +# +# plt.subplot(223) +# blh.draw(m) +# plt.title('Extended Before') +# +# m.migrad() +# +# plt.subplot(224) +# blh.draw(m) +# plt.title('Extended After') +# plt.savefig('test.png') diff --git a/experimental/rootpy_converty_hist_type.py b/experimental/rootpy_converty_hist_type.py new file mode 100644 index 00000000..263044c8 --- /dev/null +++ b/experimental/rootpy_converty_hist_type.py @@ -0,0 +1,34 @@ +''' +Created on 11 Mar 2015 + +@author: kreczko +''' + +from rootpy.io import File +from rootpy.plotting import Hist +from rootpy import asrootpy + +rootpy_hist = Hist( 100, 0, 100, type = 'F' ) + + +rootpy_hist.SetName('hist') +test_file = File('test.root', 'RECREATE') +test_file.mkdir("test") +test_file.cd('test') +rootpy_hist.Write() +test_file.Write() +test_file.Close() + +read_file = File('test.root') +folder = read_file.Get('test') +hist = folder.hist +print hist.TYPE +read_file.Close() + +hist = None + +read_file = File('test.root') +hist = read_file.Get('test/hist') +hist1 = hist.empty_clone(type='D') +print hist1.TYPE +read_file.Close() diff --git a/experimental/sums.py b/experimental/sums.py new file mode 100644 index 00000000..4caf585c --- /dev/null +++ b/experimental/sums.py @@ -0,0 +1,35 @@ +''' +Created on 23 Jan 2015 + +@author: phxlk +''' +from tools.file_utilities import read_data_from_JSON + +if __name__ == '__main__': + JSON_input_file = 'data/absolute_eta_M3_angle_bl/7TeV/HT/fit_results/central/fit_results_muon_patType1CorrectedPFMet.txt' + normalisation = read_data_from_JSON(JSON_input_file) + absolute_total_value = 0 + absolute_total_error = 0 + absolute_total_corrected_error = 0 + for sample, fit_result in normalisation.iteritems(): + print 'Calculating total # events for sample "%s"' % sample + total_events = 0 + total_error = 0 + total_corrected_error = 0 + # loop over binsZ + for result in fit_result: + value = result[0] + error = result[1] + total_events += value + total_error += error + if total_error > total_events: + total_corrected_error = total_events + else: + total_corrected_error = total_error + print 'Total number of events %d += %d (%d)' %(total_events, total_corrected_error, total_error) + absolute_total_value += total_events + absolute_total_error += total_error + absolute_total_corrected_error += total_corrected_error + print '==========================================================================' + print 'Total number of events %d += %d (%d)' %(absolute_total_value, absolute_total_corrected_error, absolute_total_error) + \ No newline at end of file diff --git a/experimental/tau_value_determination.py b/experimental/tau_value_determination.py index 42bdeeff..b592e7fc 100644 --- a/experimental/tau_value_determination.py +++ b/experimental/tau_value_determination.py @@ -23,7 +23,10 @@ from tools.file_utilities import read_data_from_JSON from tools.hist_utilities import value_error_tuplelist_to_hist from tools.Unfolding import Unfolding, get_unfold_histogram_tuple -from tools.ROOT_utililities import set_root_defaults +from tools.ROOT_utils import set_root_defaults +# from examples.Bin_Centers import nbins +from config import XSectionConfig +used_k = 2 font = {'family' : 'normal', 'weight' : 'normal', @@ -38,28 +41,59 @@ def drange( start, stop, step ): r += step def get_tau_from_global_correlation( h_truth, h_measured, h_response, h_data = None ): - tau_0 = 1e-7 - tau_max = 0.2 - number_of_iterations = 10000 + global used_k + # this gives 9.97e-05 with TUnfold + tau_0 = 1 + tau_max = 1000 + number_of_iterations = int(100) + n_toy = int(1000) # tau_step = ( tau_max - tau_0 ) / number_of_iterations optimal_tau = 0 minimal_rho = 9999 - bias_scale = 0. +# bias_scale = 0. unfolding = Unfolding( h_truth, h_measured, h_response, - method = 'RooUnfoldTUnfold', - tau = tau_0 ) + method = 'RooUnfoldSvd', + tau = tau_0, + k_value = -1, ) + data = None if h_data: - unfolding.unfold( h_data ) + data = h_data else: # closure test - unfolding.unfold( h_measured ) + data = h_measured + unfolding.unfold( data ) + # get unfolding object + tau_svd_unfold = unfolding.Impl() + # get covariance matrix + cov = tau_svd_unfold.get_data_covariance_matrix(data) # cache functions and save time in the loop - Unfold = unfolding.unfoldObject.Impl().DoUnfold - GetRho = unfolding.unfoldObject.Impl().GetRhoI - + SetTau = tau_svd_unfold.SetTau + GetCovMatrix = tau_svd_unfold.GetUnfoldCovMatrix + GetRho = tau_svd_unfold.get_global_correlation + + n_bins = h_data.nbins() + print 'k to tau' + to_return = None + for k in range(2, n_bins + 1): + tau_from_k = tau_svd_unfold.kToTau(k) + SetTau( tau_from_k ) + cov_matrix = GetCovMatrix(cov, n_toy, 1) + rho = GetRho(cov_matrix, data) + if k == used_k: + to_return = (tau_from_k, rho) + print 'k =', k, ', tau = ', tau_from_k, ', rho = ', rho, '<-- currently used' + else: + print 'k =', k, ', tau = ', tau_from_k, ', rho = ', rho + #print 'used k (=%d) to tau' % used_k + tau_from_k = tau_svd_unfold.kToTau(used_k) + #SetTau( tau_from_k ) + #cov_matrix = GetCovMatrix(cov, 10, 1) + #rho_from_tau_from_k = GetRho(cov_matrix, data) + #print "tau from k", tau_from_k + #print 'rho for tau from used k', rho_from_tau_from_k # create lists tau_values = [] rho_values = [] @@ -67,8 +101,9 @@ def get_tau_from_global_correlation( h_truth, h_measured, h_response, h_data = N add_rho = rho_values.append # for current_tau in drange(tau_0, tau_max, tau_step): for current_tau in get_tau_range( tau_0, tau_max, number_of_iterations ): - Unfold( current_tau, h_data, bias_scale ) - current_rho = GetRho( TH1F() ) + SetTau( current_tau ) + cov_matrix = GetCovMatrix(cov, n_toy, 1) + current_rho = GetRho(cov_matrix, data) add_tau( current_tau ) add_rho( current_rho ) @@ -76,10 +111,12 @@ def get_tau_from_global_correlation( h_truth, h_measured, h_response, h_data = N if current_rho < minimal_rho: minimal_rho = current_rho optimal_tau = current_tau - - return optimal_tau, minimal_rho, tau_values, rho_values + del unfolding + print 'optimal tau = %.2f' % optimal_tau + return optimal_tau, minimal_rho, tau_values, rho_values, to_return -def draw_global_correlation( tau_values, rho_values, tau, rho, channel, variable ): +def draw_global_correlation( tau_values, rho_values, tau, rho, channel, variable, tau_from_k ): + global centre_of_mass plt.figure( figsize = ( 16, 16 ), dpi = 200, facecolor = 'white' ) plt.plot( tau_values, rho_values ) plt.xscale('log') @@ -89,15 +126,24 @@ def draw_global_correlation( tau_values, rho_values, tau, rho, channel, variable ax = plt.axes() ax.annotate( r"$\tau = %.3g$" % tau, xy = ( tau, rho ), xycoords = 'data', - xytext = ( 0.0010, 0.5 ), textcoords = 'data', + xytext = ( tau*0.9, rho*1.15 ), textcoords = 'data', arrowprops = dict( arrowstyle = "fancy,head_length=0.4,head_width=0.4,tail_width=0.4", connectionstyle = "arc3" ), size = 40, ) + + ax.annotate( r"$\tau_{k} = %.3g$" % tau_from_k[0], + xy = ( tau_from_k[0], tau_from_k[1] ), xycoords = 'data', + xytext = ( tau_from_k[0], tau_from_k[1]*0.9 ), textcoords = 'data', + arrowprops = dict( arrowstyle = "fancy,head_length=0.4,head_width=0.4,tail_width=0.4", + connectionstyle = "arc3" ), + size = 40, + ) + print 'Writing to ','plots/%dTeV/tau_from_global_correlation_%s_channel_%s_DATA.png' % ( centre_of_mass, channel, variable ) if use_data: - plt.savefig( 'plots/tau_from_global_correlation_%s_channel_%s_DATA.png' % ( channel, variable ) ) + plt.savefig( 'plots/%dTeV/tau_from_global_correlation_%s_channel_%s_DATA.png' % ( centre_of_mass, channel, variable ) ) else: - plt.savefig( 'plots/tau_from_global_correlation_%s_channel_%s_MC.png' % ( channel, variable ) ) + plt.savefig( 'plots/%dTeV/tau_from_global_correlation_%s_channel_%s_MC.png' % ( centre_of_mass, channel, variable ) ) def get_tau_from_L_shape( h_truth, h_measured, h_response, h_data = None ): @@ -163,9 +209,11 @@ def draw_l_shape( l_shape, best_tau, x_value, y_value, channel, variable ): plt.savefig( 'plots/tau_from_L_shape_%s_channel_%s_MC.png' % ( channel, variable ) ) def get_data_histogram( channel, variable, met_type ): - fit_result_input = '../../data/8TeV/%(variable)s/fit_results/central/fit_results_%(channel)s_%(met_type)s.txt' + fit_result_input = 'data/M3_angle_bl/13TeV/%(variable)s/fit_results/central/fit_results_%(channel)s_%(met_type)s.txt' fit_results = read_data_from_JSON( fit_result_input % {'channel': channel, 'variable': variable, 'met_type':met_type} ) fit_data = fit_results['TTJet'] + print fit_data + print bin_edges[variable] h_data = value_error_tuplelist_to_hist( fit_data, bin_edges[variable] ) return h_data @@ -186,26 +234,43 @@ def get_tau_range( tau_min, tau_max, number_of_points ): return result; if __name__ == '__main__': + used_k = 2 set_root_defaults() use_data = True - input_file_8Tev = '/storage/TopQuarkGroup/mc/8TeV/NoSkimUnfolding/v10/TTJets_MassiveBinDECAY_TuneZ2star_8TeV-madgraph-tauola/unfolding_v10_Summer12_DR53X-PU_S10_START53_V7C-v1_NoSkim/TTJets_nTuple_53X_mc_merged_001.root' + input_files = { + 7 : '/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_5th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root', + 8 : '/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_5th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root', + 13 : 'unfolding/13TeV/unfolding_TTJets_13TeV_asymmetric.root', + } met_type = 'patType1CorrectedPFMet' # ST and HT have the problem of the overflow bin in the truth/response matrix # 7 input bins and 8 output bins (includes 1 overflow bin) - variables = ['MET', 'WPT', 'MT' , 'ST', 'HT'] - centre_of_mass = 8 - ttbar_xsection = 225.2 - luminosity = 19712 - input_file = File( input_file_8Tev ) + variables = [ + # 'MET', + # 'WPT', + # 'MT' , + 'ST', + 'HT' + ] + centre_of_mass = 13 + measurement_config = XSectionConfig(centre_of_mass) + ttbar_xsection = measurement_config.ttbar_xsection + luminosity = measurement_config.luminosity + input_file = File(input_files[centre_of_mass]) - taus_from_global_correlaton = {} - taus_from_L_shape = {} - for channel in ['electron', 'muon']: - taus_from_global_correlaton[channel] = {} - taus_from_L_shape[channel] = {} - for variable in variables: - print 'Doing variable"', variable, '" in', channel, '-channel' +# taus_from_global_correlaton = {} +# taus_from_L_shape = {} + for variable in variables: +# taus_from_global_correlaton[channel] = {} +# taus_from_L_shape[channel] = {} + for channel in ['electron', 'muon']: + if channel == 'electron': + used_k = measurement_config.k_values_electron[variable] + if channel == 'muon': + used_k = measurement_config.k_values_muon[variable] + + print 'Doing variable"', variable, '" in', channel, 'channel' h_truth, h_measured, h_response, _ = get_unfold_histogram_tuple( inputfile = input_file, @@ -223,9 +288,11 @@ def get_tau_range( tau_min, tau_max, number_of_points ): h_data = deepcopy( h_measured ) - tau, rho, tau_values, rho_values = get_tau_from_global_correlation( h_truth, h_measured, h_response, h_data ) - draw_global_correlation( tau_values, rho_values, tau, rho, channel, variable ) + print 'Getting tau' + tau, rho, tau_values, rho_values, tau_from_k = get_tau_from_global_correlation( h_truth, h_measured, h_response, h_data ) + print 'Drawing' + draw_global_correlation( tau_values, rho_values, tau, rho, channel, variable, tau_from_k ) - tau, l_curve, x, y = get_tau_from_L_shape( h_truth, h_measured, h_response, h_data ) - draw_l_shape( l_curve, tau, x, y, channel, variable ) + #tau, l_curve, x, y = get_tau_from_L_shape( h_truth, h_measured, h_response, h_data ) + #draw_l_shape( l_curve, tau, x, y, channel, variable ) diff --git a/experimental/ttbar_reco.py b/experimental/ttbar_reco.py new file mode 100644 index 00000000..f1ebfd2a --- /dev/null +++ b/experimental/ttbar_reco.py @@ -0,0 +1,257 @@ +''' +From http://arxiv.org/abs/1305.1878 +''' +import numpy as np +import ROOT as r +import math +from scipy.optimize import leastsq +mT = 172.5 +mW = 80.385 +mN = 0 +# GeV : top quark mass +# GeV : W boson mass +# GeV : neutrino mass +def UnitCircle (): + ''' + Unit circle in extended representation + ''' + return np.diag ( [1 , 1, -1] ) + +def cofactor ( A, ( i, j ) ): + '''Cofactor [i,j] of 3x3 matrix A ''' + a = A[not i:2 if i == 2 else None :2 if i == 1 else 1, + not j:2 if j == 2 else None :2 if j == 1 else 1] + return ( -1 ) ** ( i + j ) * ( a[0 , 0] * a[1 , 1] - a[1 , 0] * a[0 , 1] ) + +def R( axis , angle ): + '''Rotation matrix about x(0) ,y(1) , or z(2) axis ''' + c, s = math.cos( angle ), math.sin( angle ) + R = c * np.eye ( 3 ) + for i in [-1, 0, 1]: + R[( axis - i ) % 3 , ( axis + i ) % 3] = i * s + ( 1 - i * i ) + return R + +def Derivative (): + '''Matrix to differentiate [cos(t),sin(t) ,1] ''' + return R( 2, math.pi / 2 ).dot( np.diag ( [1 , 1, 0] ) ) + +def multisqrt ( y ): + '''Valid real solutions to y=x*x ''' + return ( [] if y < 0 else + [0] if y == 0 else + ( lambda r: [-r, r] )( math.sqrt( y ) ) ) + +def factor_degenerate ( G, zero = 0 ): + '''Linear factors of degenerate quadratic polynomial ''' + if G[0 , 0] == 0 == G[1 , 1]: + return [[G[0 , 1] , 0, G[1 , 2]] , + [0, G[0 , 1] , G[0 , 2] - G[1 , 2]]] + + swapXY = abs( G[0 , 0] ) > abs( G[1 , 1] ) + Q = G[( 1 , 0 , 2 ) , ][: , ( 1 , 0 , 2 )] if swapXY else G + Q /= Q[1 , 1] + q22 = cofactor ( Q, ( 2 , 2 ) ) + + lines = [] + if -q22 <= zero: + lines = [[Q[0, 1] , Q[1 , 1] , Q[1 , 2] + s] + for s in multisqrt ( -cofactor ( Q, ( 0 , 0 ) ) )] + else: + x0 , y0 = [ cofactor ( Q , ( i , 2 ) ) / q22 for i in [0, 1]] + lines = [[m, Q[1 , 1] , -Q[1 , 1] * y0 - m * x0] + for m in [Q[0 , 1] + s + for s in multisqrt ( -q22 )]] + + return [[L[ swapXY ], L[not swapXY ], L[2]] for L in lines ] + + +def intersections_ellipse_line( ellipse , line , zero = 1e-12 ): + '''Points of intersection between ellipse and line ''' + _, V = np.linalg.eig( np.cross( line , ellipse ).T ) + sols = sorted ( [( v.real / v[2]. real , + np.dot( line , v.real ) ** 2 + + np.dot( v.real , ellipse ).dot( v.real ) ** 2 ) + for v in V.T], + key = lambda ( s, k ): k )[:2] + return [s for s, k in sols if k < zero] + +def intersections_ellipses ( A, B, returnLines = False ): + '''Points of intersection between two ellipses ''' + LA = np.linalg + if abs( LA.det( B ) ) > abs( LA.det( A ) ): A, B = B, A + e = next( e.real for e in LA. eigvals ( LA.inv( A ).dot( B ) ) + if not e.imag ) + lines = factor_degenerate ( B - e * A ) + points = sum ( [ intersections_ellipse_line ( A, L ) + for L in lines ] , [] ) + return ( points , lines ) if returnLines else points + +class nuSolutionSet ( object ): + '''Definitions for nu analytic solution , t->b,mu ,nu ''' + + def __init__ ( self , b, mu , # Lorentz Vectors + mW2 = mW ** 2 , mT2 = mT ** 2 , mN2 = mN ** 2 ): + c = r.Math.VectorUtil.CosTheta ( b, mu ) + s = math.sqrt ( 1 - c ** 2 ) + + x0p = -( mT2 - mW2 - b.M2 () ) / ( 2 * b.E() ) + x0 = -( mW2 - mu.M2 () - mN2 ) / ( 2 * mu.E() ) + + Bb , Bm = b.Beta (), mu.Beta () + + Sx = ( x0 * Bm - mu.P() * ( 1 - Bm ** 2 ) ) / Bm ** 2 + Sy = ( x0p / Bb - c * Sx ) / s + + w = ( Bm / Bb - c ) / s + w_ = ( -Bm / Bb - c ) / s + + Om2 = w ** 2 + 1 - Bm ** 2 + eps2 = ( mW2 - mN2 ) * ( 1 - Bm ** 2 ) + x1 = Sx - ( Sx + w * Sy ) / Om2 + y1 = Sy - ( Sx + w * Sy ) * w / Om2 + Z2 = x1 ** 2 * Om2 - ( Sy - w * Sx ) ** 2 - ( mW2 - x0 ** 2 - eps2 ) + Z = math.sqrt( max ( 0, Z2 ) ) + + for item in ['b', 'mu ', 'c', 's', 'x0 ', 'x0p ', + 'Sx ', 'Sy ', 'w', 'w_ ', 'x1 ', 'y1 ', + 'Z', 'Om2 ', 'eps2 ', 'mW2 ']: + setattr ( self , item , eval( item ) ) + + @property + def K( self ): + '''Extended rotation from F' to F coord . ''' + return np.array( [[ self.c, -self.s, 0, 0], + [self.s, self.c, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1]] ) + + @property + def A_mu( self ): + '''F coord. constraint on W momentum : ellipsoid ''' + B2 = self.mu.Beta () ** 2 + SxB2 = self.Sx * B2 + F = self.mW2 - self.x0 ** 2 - self.eps2 + return np.array( [[1 - B2 , 0, 0, SxB2], + [0, 1, 0, 0], + [0, 0, 1, 0], + [SxB2, 0, 0, F]] ) + + @property + def A_b(self ): + '''F coord . constraint on W momentum : ellipsoid ''' + K, B = self.K, self.b.Beta () + mW2 , x0p = self.mW2 , self.x0p + A_b_ = np.array ([[1 -B*B, 0, 0,B*x0p], + [0, 1, 0,0], + [0, 0, 1,0], + [B*x0p , 0, 0, mW2 -x0p **2]]) + return K.dot(A_b_ ).dot(K.T) + + @property + def R_T(self ): + '''Rotation from F coord . to laboratory coord. ''' + b_xyz = self.b.x(), self.b.y(), self.b.z() + R_z = R(2, -self.mu.phi ()) + R_y = R(1, 0.5* math.pi - self.mu. theta ()) + R_x = next(R(0,- math.atan2 (z,y)) + for _,y,z in (R_y.dot(R_z.dot( b_xyz )) ,)) + return R_z.T.dot(R_y.T.dot(R_x.T)) + + @property + def H_tilde (self ): + ''' Transformation of t=[c,s ,1] to p_nu: F coord. ''' + x1 , y1 , p = self.x1 , self.y1 , self.mu.P() + Z, w, Om = self.Z, self.w, math.sqrt(self.Om2) + return np.array ([[ Z/Om , 0, x1 -p], + [w*Z/Om , 0,y1], + [0, Z,0]]) + + @property + def H(self ): + ''' Transformation of t=[c,s ,1] to p_nu: lab coord. ''' + return self.R_T.dot(self.H_tilde ) + + @property + def H_perp (self ): + ''' Transformation of t=[c,s ,1] to pT_nu : lab coord. ''' + return np.vstack([ self.H[:2] , [0, 0, 1]]) + + @property + def N(self ): + '''Solution ellipse of pT_nu : lab coord . ''' + HpInv = np.linalg.inv(self.H_perp ) + return HpInv.T.dot(UnitCircle()).dot(HpInv) + +class singleNeutrinoSolution ( object ): + '''Most likely neutrino momentum for tt --> lepton +jets ''' + def __init__ (self, b, mu, # Lorentz Vectors + (metX, metY), # Momentum imbalance + sigma2 , # Mo. imbalance unc. matrix + mW2=mW **2 , mT2=mT **2): + self.solutionSet = nuSolutionSet (b, mu , mW2 , mT2) + S2 = np.vstack ([ np.vstack ([ np.linalg.inv( sigma2 ), + [0, 0]]).T, [0, 0, 0]]) + V0 = np.outer([ metX , metY , 0], [0, 0, 1]) + deltaNu = V0 - self.solutionSet .H + + self.X = np.dot( deltaNu.T, S2 ).dot( deltaNu ) + M = next(XD + XD.T + for XD in (self.X.dot( Derivative ()) ,)) + + solutions = intersections_ellipses (M, UnitCircle ()) + self.solutions = sorted (solutions , key=self.calcX2 ) + + def calcX2 (self , t): + return np.dot(t, self.X).dot(t) + + @property + def chi2(self ): + return self.calcX2 (self.solutions [0]) + + @property + def nu(self ): + '''Solution for neutrino momentum ''' + return self.solutionSet .H.dot(self.solutions [0]) + +class doubleNeutrinoSolutions ( object ): + '''Solution pairs of neutrino momenta , tt -> leptons ''' + def __init__ (self , (b, b_), (mu , mu_), # 4- vectors + (metX , metY), # ETmiss + mW2=mW **2 , mT2=mT **2): + self.solutionSets = [ nuSolutionSet (B, M, mW2 , mT2) + for B,M in zip ((b,b_),(mu ,mu_ ))] + + V0 = np.outer ([ metX , metY , 0], [0, 0, 1]) + self.S = V0 - UnitCircle () + + N, N_ = [ss.N for ss in self.solutionSets ] + n_ = self.S.T.dot(N_ ).dot(self.S) + + v = intersections_ellipses (N, n_) + v_ = [self.S.dot(sol) for sol in v] + + if not v and leastsq : + es = [ss. H_perp for ss in self.solutionSets ] + met = np.array ([ metX , metY , 1]) + + def nus(ts ): + return tuple (e.dot ([ math.cos(t), math.sin(t), 1]) + for e, t in zip(es , ts )) + + def residuals ( params ): + return sum(nus( params ), -met )[:2] + + ts ,_ = leastsq (residuals , [0, 0], + ftol =5e-5, epsfcn =0.01) + v, v_ = [[i] for i in nus(ts )] + + for k, v in {'perp ': v, 'perp_ ': v_ , 'n_ ': n_ }. items (): + setattr (self , k, v) + + @property + def nunu_s (self ): + '''Solution pairs for neutrino momenta ''' + K, K_ = [ss.H.dot(np.linalg.inv(ss. H_perp )) + for ss in self.solutionSets ] + return [(K.dot(s), K_.dot(s_ )) + for s, s_ in zip(self.perp , self.perp_ )] \ No newline at end of file diff --git a/external/RooUnfold b/external/RooUnfold index ec05705d..5ba87eb1 160000 --- a/external/RooUnfold +++ b/external/RooUnfold @@ -1 +1 @@ -Subproject commit ec05705dbcaea3793cc82cc1c3504f90ec5dacd4 +Subproject commit 5ba87eb15ae0be09778abc73fb174af1d48acaef diff --git a/external/rootpy b/external/rootpy index b6331ccb..36446284 160000 --- a/external/rootpy +++ b/external/rootpy @@ -1 +1 @@ -Subproject commit b6331ccb0b5140fac9f8fbb11acb6db289ab7a00 +Subproject commit 36446284af5c8e3e52258360f4f36507e66f7d8e diff --git a/setup.sh b/setup.sh index c4d8bcc9..c846b062 100755 --- a/setup.sh +++ b/setup.sh @@ -67,7 +67,10 @@ pip install PrettyTable &> /dev/null success $? PrettyTable echo "Installing nose unit test framework" pip install nose &> /dev/null -success $? PrettyTable +success $? nose +echo "Installing tabulate (latex table printing, etc)" +pip install tabulate &> /dev/null +success $? tabulate echo "Installing matplotlib" #Adding freetype and libpng library and include paths from CMSSW, specific to the version but should be ok for now. @@ -88,6 +91,7 @@ if [ ! -d "$base/external/lib" ]; then mkdir $base/external/lib echo "Building RooUnfold" cd $base/external/RooUnfold/ + cmake CMakeLists.txt make -j4 success $? RooUnfold #remove tmp folder diff --git a/src/cross_section_measurement/00_pick_bins.py b/src/cross_section_measurement/00_pick_bins.py index 12eaa686..bf9a8592 100644 --- a/src/cross_section_measurement/00_pick_bins.py +++ b/src/cross_section_measurement/00_pick_bins.py @@ -44,7 +44,6 @@ from rootpy.io import File from tools.Calculation import calculate_purities, calculate_stabilities from tools.hist_utilities import rebin_2d -from config.variable_binning import bin_edges as old_binning from config import XSectionConfig from optparse import OptionParser @@ -82,10 +81,6 @@ def main(): print '-' * 120 print 'The corresponding purities and stabilities are:' for info in histogram_information: -# old_hist = rebin_2d( info['hist'], old_binning[variable], old_binning[variable] ) -# old_purities = calculate_purities( old_hist.Clone() ) -# old_stabilities = calculate_stabilities( old_hist.Clone() ) -# print_console(info, old_purities, old_stabilities, print_old = True) print_latex_table(info, variable, best_binning) print '=' * 120 diff --git a/src/unfolding_tests/get_best_regularisation.py b/src/unfolding_tests/get_best_regularisation.py new file mode 100644 index 00000000..3f8ce7e5 --- /dev/null +++ b/src/unfolding_tests/get_best_regularisation.py @@ -0,0 +1,462 @@ +''' +Created on 9 Mar 2015 + +@author: kreczko + +This script creates + - global correlation plots for tau (float) & k (int) regularisation parameters + - optimal tau & k values + - d_i plots with fits for k values + +What it needs: + - a set of four histograms: + - truth: quantity at generator level + - gen_vs_reco: quantity after selection, generated vs reco + - measured: measured (reco) quantity including fakes (background) + - data +For example config files, please have a look at config/unfolding/*.json + +usage: + python get_best_regularisation.py config.json + # for full 7 + 8 TeV stuff: + python src/unfolding_tests/get_best_regularisation.py config/unfolding/*.json -c +''' +# imports +from __future__ import division +from math import log10, pow +from optparse import OptionParser +import sys +# rootpy +from rootpy.io import File +# DailyPythonScripts +from tools.file_utilities import read_data_from_JSON, make_folder_if_not_exists +from tools.Unfolding import Unfolding, get_unfold_histogram_tuple +#from src.cross_section_measurement.lib import get_unfold_histogram_tuple +from tools.ROOT_utils import set_root_defaults, get_histogram_from_file +from config import XSectionConfig +from config.variable_binning import bin_edges +from tools.hist_utilities import value_error_tuplelist_to_hist +from tools.table import PrintTable +import matplotlib.pyplot as plt +from tools.plotting import Histogram_properties +from matplotlib import rc +from config import CMS +from config.latex_labels import variables_latex +rc('font',**CMS.font) +rc( 'text', usetex = True ) + +class RegularisationSettings(): + n_toy = int( 1000 ) + n_tau_scan_points = int( 100 ) + + def __init__( self, input_values ): + self.centre_of_mass_energy = input_values['centre-of-mass energy'] + self.measurement_config = XSectionConfig( self.centre_of_mass_energy ) + self.channel = input_values['channel'] + self.variable = input_values['variable'] + self.output_folder = input_values['output_folder'] + self.output_format = input_values['output_format'] + self.truth = input_values['truth'] + self.gen_vs_reco = input_values['gen_vs_reco'] + self.measured = input_values['measured'] + self.data = input_values['data'] + + # optional + if input_values.has_key('n_tau_scan_points'): + self.n_tau_scan_points = input_values['n_tau_scan_points'] + if input_values.has_key('n_toy'): + self.n_toy = input_values['n_toy'] + + self.__set_unfolding_histograms__() + + def __set_unfolding_histograms__( self ): + # at the moment only one file is supported for the unfolding input + files = set( [self.truth['file'], + self.gen_vs_reco['file'], + self.measured['file']] + ) + if len( files ) > 1: + print "Currently not supported to have different files for truth, gen_vs_reco and measured" + sys.exit() + + input_file = files.pop() + t, m, r, _ = get_unfold_histogram_tuple( File(input_file), + self.variable, + self.channel, + centre_of_mass = self.centre_of_mass_energy, + ttbar_xsection=self.measurement_config.ttbar_xsection, + luminosity=self.measurement_config.luminosity, + ) + self.h_truth = t + self.h_response = r + self.h_measured = m + + data_file = self.data['file'] + if data_file.endswith('.root'): + self.h_data = get_histogram_from_file(self.data['histogram'], self.data['file']) + elif data_file.endswith('.json') or data_file.endswith('.txt'): + data_key = self.data['histogram'] + # assume configured bin edges + edges = bin_edges[self.variable] + json_input = read_data_from_JSON(data_file) + if data_key == "": # JSON file == histogram + self.h_data = value_error_tuplelist_to_hist(json_input, edges) + else: + self.h_data = value_error_tuplelist_to_hist(json_input[data_key], edges) + else: + print 'Unkown file extension', data_file.split('.')[-1] + + def get_histograms( self ): + return self.h_truth, self.h_response, self.h_measured, self.h_data + +def main(): + options, input_values_sets, json_input_files = parse_options() + print 'options = ', options + use_current_k_values = options.compare + results = {} + for input_values, json_file in zip( input_values_sets, json_input_files ): + print 'Processing', json_file + regularisation_settings = RegularisationSettings( input_values ) + variable = regularisation_settings.variable + channel = regularisation_settings.channel + com = regularisation_settings.centre_of_mass_energy + if not results.has_key(com): results[com] = {} + if not results[com].has_key(variable): results[com][variable] = {} + print 'Variable = %s, channel = "%s", sqrt(s) = %d' % (variable, channel, com) + + k_results = get_best_k_from_global_correlation( regularisation_settings ) + tau_results = get_best_tau_from_global_correlation( regularisation_settings ) + + results[com][variable][channel] = (k_results, tau_results) + + plot(regularisation_settings, (k_results, tau_results), + use_current_k_values) + + table(results, use_current_k_values, options.style) + +def parse_options(): + parser = OptionParser( __doc__ ) + parser.add_option( "-c", "--compare", dest = "compare", action = "store_true", + help = "Compare to current values (k vs tau)", default = False ) + parser.add_option( "-t", "--table-style", dest = "style", default = 'simple', + help = "Style for table printing: simple|latex|twiki (default = simple)" ) + ( options, args ) = parser.parse_args() + + input_values_sets = [] + json_input_files = [] + add_set = input_values_sets.append + add_json_file = json_input_files.append + for arg in args: + input_values = read_data_from_JSON( arg ) + add_set( input_values ) + add_json_file( arg ) + + return options, input_values_sets, json_input_files + +def get_best_k_from_global_correlation( regularisation_settings ): + ''' + returns optimal_k, k_values, tau_values, rho_values + - optimal_k: k-value with lowest rho + - minimal_rho: lowest rho value + - k_values: all scanned k-values + - tau_values: tau values for all scanned k-values + - rho_values: rho values for all scanned k-values + ''' + h_truth, h_response, h_measured, h_data = regularisation_settings.get_histograms() + n_toy = regularisation_settings.n_toy + # initialise variables + optimal_k = 0 + minimal_rho = 9999 + n_bins = h_data.nbins() + k_values = [] + tau_values = [] + rho_values = [] + # first calculate one set to get the matrices + # tau = 0 is equal to k = nbins + unfolding = Unfolding( h_truth, + h_measured, + h_response, + method = 'RooUnfoldSvd', + tau = 0., # no regularisation + k_value = -1, ) + unfolding.unfold( h_data ) + # get unfolding object + svd_unfold = unfolding.Impl() + # get covariance matrix + cov = svd_unfold.get_data_covariance_matrix( h_data ) + + # cache functions and save time in the loop + SetTau = svd_unfold.SetTau + GetCovMatrix = svd_unfold.GetUnfoldCovMatrix + GetRho = svd_unfold.get_global_correlation + kToTau = svd_unfold.kToTau + add_k = k_values.append + add_tau = tau_values.append + add_rho = rho_values.append + + + # now lets loop over all possible k-values + for k in range( 2, n_bins + 1 ): + tau_from_k = kToTau( k ) + SetTau( tau_from_k ) + cov_matrix = GetCovMatrix( cov, n_toy, 1 ) + rho = GetRho( cov_matrix, h_data ) + add_k( k ) + add_tau( tau_from_k ) + add_rho( rho ) + + if rho < minimal_rho: + optimal_k = k + minimal_rho = rho + + return optimal_k, minimal_rho, k_values, tau_values, rho_values + + +def get_best_tau_from_global_correlation( regularisation_settings ): + ''' + returns optimal_tau, tau_values, rho_values + - optimal_tau: k-value with lowest rho + - minimal_rho: lowest rho value + - tau_values: all scanned tau values + - rho_values: rho values for all scanned tau-values + ''' + h_truth, h_response, h_measured, h_data = regularisation_settings.get_histograms() + n_toy = regularisation_settings.n_toy + number_of_iterations = regularisation_settings.n_tau_scan_points + tau_min = 0 + tau_max = 100 + optimal_tau = 0 + minimal_rho = 9999 + tau_values = [] + rho_values = [] + # first calculate one set to get the matrices + # tau = 0 is equal to k = nbins + unfolding = Unfolding( h_truth, + h_measured, + h_response, + method = 'RooUnfoldSvd', + tau = 0., # no regularisation + k_value = -1, ) + unfolding.unfold( h_data ) + # get unfolding object + svd_unfold = unfolding.Impl() + # get covariance matrix + cov = svd_unfold.get_data_covariance_matrix( h_data ) + + # cache functions and save time in the loop + SetTau = svd_unfold.SetTau + GetCovMatrix = svd_unfold.GetUnfoldCovMatrix + GetRho = svd_unfold.get_global_correlation + add_tau = tau_values.append + add_rho = rho_values.append + + # now lets loop over all tau-values in range + for current_tau in get_tau_range(tau_min, tau_max, number_of_iterations): + SetTau( current_tau ) + cov_matrix = GetCovMatrix(cov, n_toy, 1) + current_rho = GetRho(cov_matrix, h_data) + + add_tau( current_tau ) + add_rho( current_rho ) + + if current_rho < minimal_rho: + minimal_rho = current_rho + optimal_tau = current_tau + + return optimal_tau, minimal_rho, tau_values, rho_values + +def get_tau_range( tau_min, tau_max, number_of_points ): + # Use 3 scan points minimum + if number_of_points < 3: + number_of_points = 3 + + if tau_min == 0: + tau_min = 0.1 + + # Find the scan points + # Use equidistant steps on a logarithmic scale + step = ( log10( tau_max ) - log10( tau_min ) ) / ( number_of_points - 1 ); + for i in range ( number_of_points ): + yield pow( 10., ( log10( tau_min ) + i * step ) ); + +def plot(regularisation_settings, results, use_current_k_values = False): + variable = regularisation_settings.variable + channel = regularisation_settings.channel + com = regularisation_settings.centre_of_mass_energy + output_folder = regularisation_settings.output_folder + output_format = regularisation_settings.output_format + measurement_config = XSectionConfig(com) + + name = 'reg_param_from_global_correlation_%s_channel_%s' % (channel, variable) + hp = Histogram_properties() + hp.name = name + hp.x_axis_title = r'log($\tau$)' + hp.y_axis_title = r'$\bar{\rho}(\tau)$' + hp.title = r'global correlation for $%s$, %s channel, $\sqrt{s} = %d$ TeV' + hp.title = hp.title % (variables_latex[variable], channel, com) + + k_results, tau_results = results + optimal_tau, minimal_rho, tau_values, rho_values = tau_results + optimal_k, optimal_k_rho, k_values, k_tau_values, k_rho_values = k_results + + plt.figure( figsize = ( 16, 16 ), dpi = 200, facecolor = 'white' ) + plt.plot( tau_values, rho_values ) + plt.plot( k_tau_values, k_rho_values, 'ro' ) + + plt.title(hp.title, CMS.title) + plt.xlabel( hp.x_axis_title, CMS.x_axis_title ) + plt.ylabel( hp.y_axis_title, CMS.y_axis_title ) + plt.tick_params( **CMS.axis_label_major ) + plt.tick_params( **CMS.axis_label_minor ) + + ax = plt.axes() + + current_k, closest_tau, _, _ = get_k_tau_set(measurement_config, channel, variable, results) + current_k_rho = k_rho_values[k_values.index(current_k)] + + # first best k + tau_index = k_values.index(optimal_k) + closest_tau_best_k = k_tau_values[tau_index] + print closest_tau_best_k + ax.annotate( r"$\tau = %.2g$" % optimal_tau, + xy = ( optimal_tau, minimal_rho ), xycoords = 'data', + xytext = ( optimal_tau*0.9, minimal_rho*1.15 ), textcoords = 'data', + bbox=dict(boxstyle="round4", fc="w"), + arrowprops = dict( arrowstyle = "fancy,head_length=0.4,head_width=0.4,tail_width=0.4", + connectionstyle = "arc3" ), + size = 40, + ) + + ax.annotate( r"$\tau(k_b = %d) = %.2g$" % ( optimal_k, closest_tau_best_k ), + xy = ( closest_tau_best_k, optimal_k_rho ), xycoords = 'data', + xytext = ( closest_tau_best_k * 10, optimal_k_rho ), textcoords = 'data', + bbox=dict(boxstyle="round4", fc="w"), + arrowprops = dict( arrowstyle = "<-", + connectionstyle = "arc3", lw = 3 ), + size = 40, + ) + # then current k + if use_current_k_values: + ax.annotate( r"$\tau(k_c = %d) = %.2g$" % (current_k, closest_tau), + xy = ( closest_tau, current_k_rho ), xycoords = 'data', + xytext = ( closest_tau, current_k_rho*0.9 ), textcoords = 'data', + bbox=dict(boxstyle="round4", fc="w"), + arrowprops = dict( arrowstyle = "<-", + connectionstyle = "arc3", lw = 3 ), + size = 40, + ) + + plt.xscale('log') + make_folder_if_not_exists(output_folder) + for f in output_format: + plt.savefig( output_folder + '/' + hp.name + '.' + f ) + + + +def table(result_dict, use_current_k_values = False, style = 'simple'): + ''' + result_dict has the form + { + centre-of-mass energy : + { + variable : + { + channel : (k_results, tau_results) + } + } + } + _results are of the form: + best_, best__rho, s, rhos + ''' + for com in result_dict.keys(): + # step 1: group the results by centre of mass energy + headers = [] + if use_current_k_values: + headers = ['Variable', 'current k', 'closest tau', 'best tau', 'best k'] + else: + headers = ['Variable', 'best k', 'rho (best k)', 'best tau', 'rho (best tau)'] + data = [] + measurement_config = XSectionConfig(com) + for variable in result_dict[com].keys(): + has_both_channels = len(result_dict[com][variable]) == 2 + # step 2: if have electron and muon channel, group them: electron (muon) + if has_both_channels: + electron_results = result_dict[com][variable]['electron'] + muon_results = result_dict[com][variable]['muon'] + + entry = [] + if use_current_k_values: + electron_set = get_k_tau_set(measurement_config, 'electron', + variable, electron_results) + muon_set = get_k_tau_set(measurement_config, 'muon', + variable, muon_results) + + entry = [variable, + '%d (%d)' % (electron_set[0], muon_set[0]), + '%.1f (%.1f)' % (electron_set[1], muon_set[1]), + '%.1f (%.1f)' % (electron_set[2], muon_set[2]), + '%d (%d)' % (electron_set[3], muon_set[3]), + ] + else: + entry = [variable, + '%d (%d)' % (electron_results[0][0], muon_results[0][0]), + '%.1f (%.1f)' % (electron_results[0][1], muon_results[0][1]), + '%.1f (%.1f)' % (electron_results[1][0], muon_results[1][0]), + '%.1f (%.1f)' % (electron_results[1][1], muon_results[1][1]), + ] + + data.append(entry) + else: + channel = result_dict[com][variable].keys()[0] + results = result_dict[com][variable][channel] + + if use_current_k_values: + result_set = get_k_tau_set(measurement_config, channel, + variable, results) + entry = [variable, + '%d' % result_set[0], + '%.1f' % result_set[1], + '%.1f' % result_set[2], + '%d' % result_set[3], + ] + else: + entry = [variable, + '%d' % results[0][0], + '%.1f' % results[0][1], + '%.1f' % results[1][0], + '%.1f' % results[1][1], + ] + + data.append(entry) + + table = PrintTable(data, headers) + + print 'Printing table for sqrt(s) = %d TeV' % com + if style == 'simple': + print table.simple() + elif style == 'latex': + print table.latex() + elif style == 'twiki': + print table.twiki() + else: + print 'unknown printing style "%s"' % style + +def get_k_tau_set(config, channel, variable, results): + k_results, tau_results = results + optimal_tau, _, tau_values, rho_values = tau_results + optimal_k, _, k_values, _, k_rho_values = k_results + + current_k = config.k_values_electron[variable] + rho_index = k_values.index(current_k) + current_k_rho = k_rho_values[rho_index] + + tau_index = min(range(len(rho_values)), + key = lambda i: abs(rho_values[i] - current_k_rho)) + closest_tau = tau_values[tau_index] + + return current_k, closest_tau, optimal_tau, optimal_k + + +if __name__ == '__main__': + set_root_defaults() + main() diff --git a/test/tools_Calculation.py b/test/tools_Calculation.py index 6165f3bb..2415b6f2 100644 --- a/test/tools_Calculation.py +++ b/test/tools_Calculation.py @@ -76,7 +76,8 @@ def test_random_elipse_purity( self ): # purities should always be above 0 and below ~0.5 for p in purities: self.assertGreater( p, 0 ) - self.assertLess( p, 0.5 ) + # allow for 10% error due to randomness + self.assertLess( p, 0.5 + 0.05) def test_random_elipse_stability( self ): stabilities = calculate_stabilities( self.random_elipse ) @@ -85,7 +86,8 @@ def test_random_elipse_stability( self ): # stabilities should always be above 0 and below ~0.6 for s in stabilities: self.assertGreater( s, 0 ) - self.assertLess( s, 0.6 ) + # allow for 10% error due to randomness + self.assertLess( s, 0.6 + 0.06) def test_pre_calculated_purity( self ): purities = calculate_purities( self.pre_calculated ) diff --git a/test/tools_Unfolding.py b/test/tools_Unfolding.py index 52f2b022..2ad54647 100644 --- a/test/tools_Unfolding.py +++ b/test/tools_Unfolding.py @@ -40,31 +40,59 @@ def setUp( self ): k_value = self.k_value, method = self.unfold_method ) + + tau_unfolding_object = Unfolding( h_truth, + h_measured, + h_response, + tau=100, + k_value= -1, + method='RooUnfoldSvd') self.dict[channel][variable] = {'h_truth' : h_truth, 'h_measured' : h_measured, 'h_response' : h_response, - 'unfolding_object' : unfolding_object + 'unfolding_object' : unfolding_object, + 'tau_unfolding_object': tau_unfolding_object, } def tearDown( self ): pass - def test_closure( self ): +# def test_closure( self ): +# for channel in self.channels: +# for variable in self.variables: +# # closure test +# unfolded_result = hist_to_value_error_tuplelist( self.dict[channel][variable]['unfolding_object'].closureTest() ) +# truth = hist_to_value_error_tuplelist( self.dict[channel][variable]['h_truth'] ) +# # the difference between the truth and unfolded result should be within the unfolding error +# for (value, error), (true_value, _) in zip(unfolded_result, truth): +# self.assertAlmostEquals(value, true_value, delta = error) + + def test_invalid_zero_data( self ): + variable = 'MET' + channel = 'electron' + pseudo_data = value_error_tuplelist_to_hist( [(0,0)]*( len( bin_edges[variable] ) - 1 ), bin_edges[variable] ) + self.assertRaises(ValueError, self.dict[channel][variable]['unfolding_object'].unfold, (pseudo_data)) + + def test_tau_closure(self): for channel in self.channels: for variable in self.variables: - # closure test - unfolded_result = hist_to_value_error_tuplelist( self.dict[channel][variable]['unfolding_object'].closureTest() ) + data = self.dict[channel][variable]['h_measured'] truth = hist_to_value_error_tuplelist( self.dict[channel][variable]['h_truth'] ) + unfolded_result = hist_to_value_error_tuplelist(self.dict[channel][variable]['tau_unfolding_object'].unfold(data)) # the difference between the truth and unfolded result should be within the unfolding error for (value, error), (true_value, _) in zip(unfolded_result, truth): self.assertAlmostEquals(value, true_value, delta = error) + print value, '+-', error, ' true:', true_value - def test_invalid_zero_data( self ): - variable = 'MET' - channel = 'electron' - pseudo_data = value_error_tuplelist_to_hist( [(0,0)]*( len( bin_edges[variable] ) - 1 ), bin_edges[variable] ) - self.assertRaises(ValueError, self.dict[channel][variable]['unfolding_object'].unfold, (pseudo_data)) + def test_k_to_tau(self): + data = self.dict['electron']['MET']['h_measured'] + tau_unfolding_object = self.dict['electron']['MET']['tau_unfolding_object'] + # first we need to unfold to get the matrix + tau_unfolding_object.unfold(data) + # next we need to get the actual RooUnfold object + tau = tau_unfolding_object.Impl().kToTau(self.k_value) + self.assertAlmostEqual(tau, 19., delta = 1) if __name__ == "__main__": set_root_defaults() diff --git a/tools/ROOT_utils.py b/tools/ROOT_utils.py index b79d9f14..f6b5ccfd 100644 --- a/tools/ROOT_utils.py +++ b/tools/ROOT_utils.py @@ -5,8 +5,7 @@ ''' from rootpy.logger import logging from rootpy.io import File -from rootpy.plotting import Hist -from ROOT import gROOT +from ROOT import gROOT, TH1F gcd = gROOT.cd from config.summations_common import b_tag_bins_inclusive, b_tag_summations from config.summations_common import b_tag_bins_exclusive @@ -37,7 +36,13 @@ def get_histogram_from_file( histogram_path, input_file ): root_histogram.Add( hist ) gcd() - histogram = root_histogram.Clone() + histogram = None + # change from float to double + if root_histogram.TYPE == 'F': + histogram = root_histogram.empty_clone(type='D') + histogram.Add(root_histogram) + else: + histogram = root_histogram.Clone() root_file.Close() return histogram @@ -142,6 +147,8 @@ def set_root_defaults( set_batch = True, msg_ignore_level = 1001 ): gROOT.SetBatch( set_batch ) # ignore warnings gROOT.ProcessLine( 'gErrorIgnoreLevel = %d;' % msg_ignore_level ) + # turn of the stupid ROOT pointer handling, seriously + TH1F.AddDirectory(False) def root_mkdir(file_handle, path): ''' diff --git a/tools/Unfolding.py b/tools/Unfolding.py index d6e922cb..791f11aa 100644 --- a/tools/Unfolding.py +++ b/tools/Unfolding.py @@ -41,14 +41,12 @@ def __init__( self, self.fakes = fakes self.response = response self.data = None - self.unfolded_closure = None self.unfolded_data = None self.unfoldObject = None self.unfoldResponse = None - self.closure_test = None self.verbose = verbose - self.tau = tau - self.k_value = k_value + self.tau = float(tau) + self.k_value = int(k_value) self.n_toy = n_toy self.Bayes_n_repeat = Bayes_n_repeat self.Hreco = Hreco @@ -66,7 +64,11 @@ def setup_unfolding ( self, data ): elif self.method == 'RooUnfoldInvert': self.unfoldObject = RooUnfoldInvert ( self.unfoldResponse, self.data ) elif self.method == 'RooUnfoldSvd': - self.unfoldObject = RooUnfoldSvd( self.unfoldResponse, self.data, self.k_value, self.n_toy ) + if self.k_value > 0: + self.unfoldObject = RooUnfoldSvd( self.unfoldResponse, self.data, self.k_value, self.n_toy ) + else: + if self.tau >= 0: + self.unfoldObject = RooUnfoldSvd( self.unfoldResponse, self.data, self.tau, self.n_toy ) elif self.method == 'RooUnfoldTUnfold': self.unfoldObject = RooUnfoldTUnfold ( self.unfoldResponse, data ) if self.tau >= 0: @@ -146,6 +148,8 @@ def test_regularisation ( self, data, k_max ): raise ValueError( 'Unfolding method "%s" is not supported for regularisation parameter tests. Please use RooUnfoldSvd.' % ( self.method ) ) def unfold( self, data ): + if data is None: + raise ValueError('Data histogram is None') have_zeros = [value == 0 for value,_ in hist_to_value_error_tuplelist( data )] if not False in have_zeros: raise ValueError('Data histograms contains only zeros') @@ -159,40 +163,7 @@ def unfold( self, data ): return self.unfolded_data def closureTest( self ): - if not self.closure_test: - if not self.unfoldResponse: - self.unfoldResponse = self._makeUnfoldResponse() - if self.method == 'RooUnfoldBayes': - self.closure_test = RooUnfoldBayes ( self.unfoldResponse, self.measured, self.Bayes_n_repeat ) - elif self.method == 'RooUnfoldBinByBin': - self.closure_test = RooUnfoldBinByBin ( self.unfoldResponse, self.measured ) - elif self.method == 'RooUnfoldInvert': - self.closure_test = RooUnfoldInvert ( self.unfoldResponse, self.measured ) - elif self.method == 'RooUnfoldSvd': - self.closure_test = RooUnfoldSvd( self.unfoldResponse, self.measured, self.k_value, self.n_toy ) - elif self.method == 'RooUnfoldTUnfold': - self.closure_test = RooUnfoldTUnfold ( self.unfoldResponse, self.measured ) - if self.tau >= 0: - self.closure_test.FixTau( self.tau ) - elif self.method == 'TSVDUnfold': - new_measured = Hist( list( self.measured.xedges() ), type = 'D' ) - new_measured.Add( self.measured ) - new_truth = Hist( list( self.truth.xedges() ), type = 'D' ) - new_truth.Add( self.truth ) - - if self.fakes: - new_fakes = Hist( list ( self.fakes.xedges() ), type = 'D' ) - new_fakes.Add ( self.fakes ) - new_measured = new_measured - new_fakes - - new_response = Hist2D( list( self.response.xedges() ), list( self.response.yedges() ), type = 'D' ) - new_response.Add( self.response ) - self.closure_test = TSVDUnfold( new_measured, new_measured, new_truth, new_response ) - if self.method == 'TSVDUnfold': - self.unfolded_closure = asrootpy( self.closure_test.Unfold( self.k_value ) ) - else: - self.unfolded_closure = asrootpy( self.closure_test.Hreco( self.Hreco ) ) - return self.unfolded_closure + return self.unfold(self.measured) def _makeUnfoldResponse( self ): if self.fakes: @@ -269,6 +240,9 @@ def get_unfolded_data_errors( self ): decomposition.Solve( unfolded_errors ); return unfolded_errors + + def Impl(self): + return self.unfoldObject.Impl() def get_unfold_histogram_tuple( inputfile, @@ -279,7 +253,7 @@ def get_unfold_histogram_tuple( ttbar_xsection = 245.8, luminosity = 19712, load_fakes = False, - scale_to_lumi = True, + scale_to_lumi = False, ): folder = None h_truth = None @@ -294,7 +268,6 @@ def get_unfold_histogram_tuple( h_truth = asrootpy( folder.truth.Clone() ) h_measured = asrootpy( folder.measured.Clone() ) - # response matrix is always without fakes # fake subtraction from measured is performed automatically in RooUnfoldSvd (h_measured - h_response->ProjectionX()) # or manually for TSVDUnfold diff --git a/tools/table.py b/tools/table.py new file mode 100644 index 00000000..d4ac7600 --- /dev/null +++ b/tools/table.py @@ -0,0 +1,72 @@ +''' +Created on 25 Mar 2015 + +@author: kreczko +''' +from tabulate import tabulate + +class PrintTable(): + ''' + Class for the creation and printing of tables + Example usage: + data = [ + ['DailyPythonScripts', 'Python', 'last parts of analysis'], + ['NTupleProduction', 'C++/CMSSW', 'miniAOD -> ntuple'], + ] + headers = ['software package', 'main language', 'comment'] + t = PrintTable(data, headers) + # simple nice table (for console) + print t.simple() + # table to put into a twiki + print t.twiki() + # table to put into a latex document + print t.latex() + ''' + + + def __init__(self, data, headers): + ''' + @param data: a list of columns + @param headers: a list of headers (must be of column length) + ''' + self.data = data + self.headers = headers + + def twiki(self): + ''' + return the table as a string in twiki format + ''' + # add bold for twiki headers + headers = ['*' + h + '*' for h in self.headers] + # use Emacs org-mode as baseline + table = tabulate(self.data, headers = headers, tablefmt='orgtbl') + table = str(table).split('\n') + twiki_table = [] + for line in table: + # remove header delimiter + if not line.startswith('|--'): + twiki_table.append(line) + + return ''.join([line + '\n' for line in twiki_table]) + + + def simple(self): + ''' + return the table as a pretty string + ''' + table = tabulate(self.data, headers = self.headers, tablefmt='simple') + return str(table) + + def latex(self): + ''' + return the table as a string in latex format + ''' + table = tabulate(self.data, headers = self.headers, tablefmt='latex') + return str(table) + + +if __name__ == '__main__': + data = [['MET','X', ' '], ['MT',' ', 'x']] + headers = ['variable', 'fun', 'not fun'] + table = PrintTable(data, headers) + print table.twiki()