Plotting Tutorial

CMakeLists.txt

@@ -115,6 +115,7 @@ install(FILES "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.MaCh3Tutorial
 
 ############################  Install  ####################################
 install(DIRECTORY Inputs DESTINATION ${CMAKE_BINARY_DIR})
+install(DIRECTORY plotting DESTINATION ${CMAKE_BINARY_DIR})
 install(DIRECTORY Scripts DESTINATION ${CMAKE_BINARY_DIR})
 
 include(CMakePackageConfigHelpers)

README.md

@@ -25,6 +25,13 @@ docker pull ghcr.io/mach3-software/mach3tutorial:alma9latest
 ```
 To reed more how to use containers check our wiki [here](https://github.com/mach3-software/MaCh3/wiki/12.-Containers)
 
+## How to run LLH scan
+To run an LLH scan simply do
+```
+./bin/MCMCTutorial Inputs/ManagerTest.yaml
+```
+These are very useful for validations and to compare with other fitters
+
 ## How to run MCMC
 To run MCMC simply
 ```
@@ -34,10 +41,80 @@ Congratulations! 🎉
 You have just completed finished you first MCMC chain.
 
 ## How to Plot?
+
+There are a number of apps included to make plots from the results of your fits, llh scans etc. You can find more details on them and how they work in the main MaCh3 wiki [here](https://github.com/mach3-software/MaCh3/wiki). There you will also find some instructions on how you can write yor own plotting scripts. 
+
+The plotting library is configured using yaml files. You can see some examples of such config files in the plotting directory, and a detailed explanation of them is given in [the wiki](https://github.com/mach3-software/MaCh3/wiki).
+
+Some examples on how to make some "standard" plots are given below.
+
+### MCMC Chain
 You can produce simple plots with
 ```
 ./bin/ProcessMCMC bin/TutorialDiagConfig.yaml Test.root
 ```
+where Test.root is the output of running MCMCTutorial as described [here](#how-to-run-mcmc)
+
+You can then take the output of running ProcessMCMC which will be called something like <inputName>_Process.root, and make fancier error plots from it using the `GetPostfitParamPlots` app like 
+
+```
+GetPostfitParamPlots Test_Process.root
+```
+
+### LLH Scans
+
+You can plot the results of an LLH scan using the aptly named PlotLLH app like so
+
+```
+PlotLLH LLH_Test.root
+```
+
+where LLH_Test.root is the result of running the LLH scan as described [here](#how-to-run-llh-scan).
+
+### Plotting with Python
+
+If you have installed the python interface for MaCh3 as described [here](https://github.com/mach3-software/MaCh3?tab=readme-ov-file#python) then you can also use the provided python plotting module. The details on how to write custom python scripts using this plotting module are detailed in [the wiki](https://github.com/mach3-software/MaCh3/wiki/15.-Plotting#custom-plotting-scripts). Here we will walk you through some example scripts.
+
+For the examples here, we will use matplotlib and numpy. These can be installed using the provided [requirements.txt](requirements.txt) by doing 
+
+```
+pip install -r requirements.txt
+```
+
+but note that these are just what is used for an example for the purpose of this tutorial. When making your own plots you are totally free to use whatever plotting libraries you like!
+
+You can use the example script [MCMCPlotting-tutorial.py](plotting/MCMC-plotting-tutorial.py) to plot the raw MCMC chain values that were obtained in the [how to run MCMC](#how-to-run-mcmc) section above by doing 
+
+```
+python plotting/MCMC-plotting-tutorial.py Test.root
+```
+
+This will give you some plots that look something like 
+
+<img width="350" alt="MCMC example" src="https://github.com/user-attachments/assets/bdb1792f-3d52-4ea3-8eb8-0e6cf7b9119a">
+
+After you have run this chain through the MCMC processor as described in the [How To Plot?](#how-to-plot) section, you can pass the processed file to [1DPosterior-tutorial.py](plotting/1DPosterior-tutorial.py) like
+
+```
+python plotting/1DPosterior-tutorial.py Test_Process.root
+```
+
+which will plot the projected one dimensional posteriors which should look something like 
+
+<img width="350" alt="1dPosterior example" src="https://github.com/user-attachments/assets/4ddf3abb-9794-4f21-ae9b-862718c2ff57">
+
+
+Similarly you can use [LLHScan-plotting-tutorial.py](plotting/LLHScan-plotting-tutorial.py) to plot the LLH scans you made in the [How to run LLH scan](#how-to-run-llh-scan) section like
+
+```
+python plotting/LLHScan-plotting-tutorial.py LLH_Test.root
+```
+
+which will give you some plots that look something like 
+
+<img width="350" alt="LLH scan example" src="https://github.com/user-attachments/assets/f16ad571-68da-42e3-ae6b-e984d03a58c3">
+
+
 
 ## How to Setup your analysis
 In the next step you gonna modify analysis setup and repeat steps.

Tutorial/CMakeLists.txt

@@ -2,6 +2,7 @@ add_custom_target(TutorialApps)
 
 foreach(app
     MCMCTutorial
+    LLHScanTutorial
     )
     add_executable(${app} ${app}.cpp)
     target_link_libraries(${app} MaCh3Validations::ValidationsUtils)

Tutorial/LLHScanTutorial.cpp

@@ -0,0 +1,31 @@
+// MaCh3 spline includes
+#include "mcmc/MaCh3Factory.h"
+#include "covariance/covarianceXsec.h"
+
+int main(int argc, char *argv[]){
+
+  MaCh3Utils::MaCh3Usage(argc, argv);
+  // Initialise manger responsible for config handling
+  manager *FitManager = new manager(argv[1]);
+
+  std::vector<std::string> xsecCovMatrixFile = FitManager->raw()["General"]["Systematics"]["XsecCovFile"].as<std::vector<std::string>>();
+  // Initialise covariance class reasonable for Systematics
+  covarianceXsec* xsec = new covarianceXsec(xsecCovMatrixFile, "xsec_cov");
+
+  std::vector<std::string> OscMatrixFile = FitManager->raw()["General"]["Systematics"]["OscCovFile"].as<std::vector<std::string>>();
+  covarianceOsc* osc = new covarianceOsc(OscMatrixFile, "osc_cov");
+
+  // Create MCMC Class
+  std::unique_ptr<FitterBase> MaCh3Fitter = std::make_unique<mcmc>(FitManager);
+  // Add covariance to MCM
+  MaCh3Fitter->addSystObj(xsec);
+  MaCh3Fitter->addSystObj(osc);
+
+  // Run MCMCM
+  MaCh3Fitter->RunLLHScan();
+
+  delete FitManager;
+  delete xsec;
+  MaCh3Fitter.reset();
+  return 0;
+}

plotting/1DPosterior-tutorial.py

@@ -0,0 +1,43 @@
+## import the plotting manager from the plotting library
+from pyMaCh3.pyMaCh3.plotting import PlottingManager
+
+## this is needed in order to get the command line arguments
+import sys
+
+## get the matplotlib stuff
+from matplotlib import pyplot as plt
+import matplotlib.collections as mcoll
+from matplotlib import cm
+import matplotlib.backends.backend_pdf
+
+## some other python libraries
+import numpy as np
+import math as m
+import imageio
+
+def plot_1d_posteriors(man):
+    pdf = matplotlib.backends.backend_pdf.PdfPages("1d-posteriors.pdf")
+
+    # loop through all parameters the manager knows about
+    for i, param in enumerate( man.input().get_known_parameters() ):
+
+        fig = plt.figure()
+        plt.title(man.style().prettify_parameter_name(param))
+
+        # go through all files specified in cmd line
+        for file_id in range(man.input().get_n_input_files()):
+            posterior = man.input().get_1d_posterior(file_id, param)
+            plt.plot(posterior[0], posterior[1], label = man.get_file_label(file_id))
+
+        plt.ylabel("N Steps")
+        plt.legend()
+
+        pdf.savefig( fig )
+
+    pdf.close()
+
+# make a PlottingManager and initialise using the command line arguments
+manager = PlottingManager()
+manager.parse_inputs(sys.argv)
+
+plot_1d_posteriors(manager)

plotting/LLHScan-plotting-tutorial.py

@@ -0,0 +1,42 @@
+## import the plotting manager from the plotting library
+from pyMaCh3.pyMaCh3.plotting import PlottingManager
+
+## this is needed in order to get the command line arguments
+import sys
+
+## get the matplotlib stuff
+from matplotlib import pyplot as plt
+import matplotlib.backends.backend_pdf
+
+## some other python libraries
+import numpy as np
+import math as m
+
+def test_LLH_scans(man):
+    ''' Function to plot log likelihood scans '''
+
+    pdf = matplotlib.backends.backend_pdf.PdfPages("LLH-scans.pdf")
+
+    # loop through all parameters the manager knows about
+    for i, param in enumerate( man.input().get_known_parameters() ):
+
+        fig = plt.figure()
+        plt.title(man.style().prettify_parameter_name(param))
+
+        # go through all files specified in cmd line
+        for file_id in range(man.input().get_n_input_files()):
+            llh_scan = man.input().get_llh_scan(file_id, param, "total")
+            plt.plot(llh_scan[0], llh_scan[1], label = man.get_file_label(file_id))
+
+        plt.ylabel("LLH")
+        plt.legend()
+
+        pdf.savefig( fig )
+
+    pdf.close()
+
+# make a PlottingManager and initialise using the command line arguments
+manager = PlottingManager()
+manager.parse_inputs(sys.argv)
+
+test_LLH_scans(manager)

plotting/MCMC-plotting-tutorial.py

@@ -0,0 +1,106 @@
+## import the plotting manager from the plotting library
+from pyMaCh3.pyMaCh3.plotting import PlottingManager
+
+## this is needed in order to get the command line arguments
+import sys
+
+## get the matplotlib stuff
+from matplotlib import pyplot as plt
+import matplotlib.collections as mcoll
+from matplotlib import cm
+import matplotlib.backends.backend_pdf
+
+## some other python libraries
+import numpy as np
+import math as m
+import imageio
+
+cmap = cm.inferno
+
+def make_segments(x, y):
+    """
+    Create list of line segments from x and y coordinates, in the correct format
+    for LineCollection: an array of the form numlines x (points per line) x 2 (x
+    and y) array
+    """
+
+    points = np.array([x, y]).T.reshape(-1, 1, 2)
+    segments = np.concatenate([points[:-1], points[1:]], axis=1)
+    return segments
+
+def test_MCMC(man):
+
+    ## make the plot
+    fig, ax = plt.subplots() 
+
+    ## The number of steps in the chain
+    n_steps = man.input().get_n_MCMC_entries(0)
+
+    pdf = matplotlib.backends.backend_pdf.PdfPages("MCMC-plots.pdf")
+
+    for param_1 in man.input().get_tagged_parameters(["oscillation"], "any"):
+
+        ## check that param_1 actually has MCMC steps
+        if not param_1 in man.input().get_known_MCMC_parameters(0): 
+            continue
+
+        for param_2 in man.input().get_tagged_parameters(["oscillation"], "any"):
+
+            ## check that param_2 actually has MCMC steps
+            if not param_2 in man.input().get_known_MCMC_parameters(0): 
+                continue
+
+            if param_1 == param_2: 
+                continue
+
+            print ( "on parameters {} and {}".format(param_1, param_2) )
+
+            ## hold the values of each of the 2 parameters
+            ax1_vals = []
+            ax2_vals = []
+
+            for step in range(0, n_steps):
+                if ( step % m.floor(n_steps / 5) == 0 ):
+                    print ( "Step {} / {} :: {} %".format(step, n_steps, 100.0 * float(step) / float(n_steps)) )
+
+                man.input().get_MCMC_entry(0, step)
+
+                ax1_vals.append( man.input().get_MCMC_value(0, param_1))
+                ax2_vals.append( man.input().get_MCMC_value(0, param_2))
+
+
+            ## make line segments and colour gradient values
+            segments = make_segments( ax1_vals, ax2_vals )
+            colours = np.arange( n_steps )
+
+            ## make the collection of line segments
+            lc = mcoll.LineCollection(segments, linewidth=1, alpha=0.75, array=colours, cmap=cmap)
+
+
+            plt.cla()
+            ax.add_collection(lc)
+            ax.set_facecolor('k')
+
+            plt.xlabel(man.style().prettify_parameter_name(param_1))
+            plt.ylabel(man.style().prettify_parameter_name(param_2))
+
+            ax1_array = np.array(ax1_vals)
+            ax2_array = np.array(ax2_vals)
+
+            plt.xlim(0.95 * ax1_array.min(), 1.05 * ax1_array.max())
+            plt.ylim(0.95 * ax2_array.min(), 1.05 * ax2_array.max())
+
+            cbar = plt.colorbar(lc, label="Step")
+            plt.plot()
+
+            pdf.savefig( fig )
+
+            cbar.remove()
+
+    pdf.close()
+
+# make a PlottingManager and initialise
+manager = PlottingManager()
+manager.parse_inputs(sys.argv)
+
+test_MCMC(manager)