#!/usr/bin/env python3 # Copyright (c) 2019, Anthony Latorre # # This program is free software: you can redistribute it and/or modify it # under the terms of the GNU General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) # any later version. # # This program is distributed in the hope that it will be useful, but WITHOUT # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or # FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for # more details. # # You should have received a copy of the GNU General Public License along with # this program. If not, see . """ Script to plot reconstructed quantities for instrumentals. To run it just run: $ ./plot-dc [list of fit results] This will produce corner plots showing the distribution of the high level variables used in the contamination analysis for all the different instrumental backgrounds and external muons. """ from __future__ import print_function, division import numpy as np from scipy.stats import iqr, poisson from matplotlib.lines import Line2D from scipy.stats import iqr, norm, beta from scipy.special import spence from itertools import chain particle_id = {20: 'e', 22: r'\mu'} def plot_hist2(df, muons=False): for id, df_id in sorted(df.groupby('id')): if id == 20: plt.subplot(2,3,1) elif id == 22: plt.subplot(2,3,2) elif id == 2020: plt.subplot(2,3,4) elif id == 2022: plt.subplot(2,3,5) elif id == 2222: plt.subplot(2,3,6) if muons: plt.hist(np.log10(df_id.ke.values/1000), bins=np.linspace(0,4.5,100), histtype='step') plt.xlabel("log10(Energy (GeV))") else: bins = np.logspace(np.log10(20),np.log10(10e3),21) plt.hist(df_id.ke.values, bins=bins, histtype='step') plt.gca().set_xscale("log") major = np.array([10,100,1000,10000]) minor = np.unique(list(chain(*list(range(i,i*10,i) for i in major[:-1])))) minor = np.setdiff1d(minor,major) plt.gca().set_xticks(major) plt.gca().set_xticks(minor,minor=True) plt.xlabel("Energy (MeV)") plt.title('$' + ''.join([particle_id[int(''.join(x))] for x in grouper(str(id),2)]) + '$') if len(df): plt.tight_layout() def plot_hist(df, muons=False): for id, df_id in sorted(df.groupby('id')): if id == 20: plt.subplot(3,4,1) elif id == 22: plt.subplot(3,4,2) elif id == 2020: plt.subplot(3,4,5) elif id == 2022: plt.subplot(3,4,6) elif id == 2222: plt.subplot(3,4,7) elif id == 202020: plt.subplot(3,4,9) elif id == 202022: plt.subplot(3,4,10) elif id == 202222: plt.subplot(3,4,11) elif id == 222222: plt.subplot(3,4,12) if muons: plt.hist(np.log10(df_id.ke.values/1000), bins=np.linspace(0,4.5,100), histtype='step') plt.xlabel("log10(Energy (GeV))") else: plt.hist(df_id.ke.values, bins=np.linspace(20,10e3,100), histtype='step') plt.xlabel("Energy (MeV)") plt.title(str(id)) if len(df): plt.tight_layout() if __name__ == '__main__': import argparse import numpy as np import pandas as pd import sys import h5py from sddm.plot_energy import * from sddm.plot import * from sddm import setup_matplotlib parser = argparse.ArgumentParser("plot fit results") parser.add_argument("filenames", nargs='+', help="input files") parser.add_argument("--save", action='store_true', default=False, help="save corner plots for backgrounds") args = parser.parse_args() setup_matplotlib(args.save) import matplotlib.pyplot as plt # Loop over runs to prevent using too much memory evs = [] data_filenames = [filename for filename in args.filenames if 'SNO' in filename] mc_filenames = [filename for filename in args.filenames if 'SNO' not in filename] if len(data_filenames): rhdr = pd.concat([read_hdf(filename, "rhdr").assign(filename=filename) for filename in data_filenames],ignore_index=True) for run, df in rhdr.groupby('run'): evs.append(get_events(df.filename.values, merge_fits=True)) if len(mc_filenames): for filename in mc_filenames: evs.append(get_events([filename], merge_fits=True, mc=True)) ev = pd.concat([ev for ev in evs if len(ev) > 0]) ev = ev[ev.prompt & ~np.isnan(ev.fmin)] ev = ev[ev.ke > 20] #with pd.option_context('display.max_rows', None, 'display.max_columns', None): # print("Noise events") # print(ev[ev.noise][['psi','x','y','z','id1','id2']]) # print("Muons") # print(ev[ev.muon][['psi','r','id1','id2','id3','energy1','energy2','energy3']]) # print("Neck") # print(ev[ev.neck & ev.psi < 6][['psi','r','id1','cos_theta']]) # print("Flashers") # print(ev[ev.flasher & ev.udotr > 0]) # print("Signal") # print(ev[ev.signal]) # save as PDF b/c EPS doesn't support alpha values if args.save: plot_corner_plot(ev[ev.breakdown],"Breakdowns",save="breakdown_corner_plot") plot_corner_plot(ev[ev.muon],"Muons",save="muon_corner_plot") plot_corner_plot(ev[ev.flasher],"Flashers",save="flashers_corner_plot") plot_corner_plot(ev[ev.neck],"Neck",save="neck_corner_plot") plot_corner_plot(ev[ev.noise],"Noise",save="noise_corner_plot") plot_corner_plot(ev[ev.signal],"Signal",save="signal_corner_plot") else: plot_corner_plot(ev[ev.breakdown],"Breakdowns") plot_corner_plot(ev[ev.muon],"Muons") plot_corner_plot(ev[ev.flasher],"Flashers") plot_corner_plot(ev[ev.neck],"Neck") plot_corner_plot(ev[ev.noise],"Noise") plot_corner_plot(ev[ev.signal],"Signal") fig = plt.figure() plot_hist2(ev[ev.flasher]) despine(fig,trim=True) plt.suptitle("Flashers") fig = plt.figure() plot_hist2(ev[ev.muon],muons=True) despine(fig,trim=True) plt.suptitle("Muons") plt.show() ='n120' href='#n120'>120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444
/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
 *
 * This program is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option)
 * any later version.

 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details.

 * You should have received a copy of the GNU General Public License along with
 * this program. If not, see <https://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <math.h>
#include "optics.h"
#include "quantum_efficiency.h"
#include "solid_angle.h"
#include "pdg.h"
#include "vector.h"
#include "electron.h"
#include "sno.h"
#include "scattering.h"
#include "pmt_response.h"
#include "misc.h"
#include <gsl/gsl_integration.h>
#include <math.h> /* for fmax() */
#include <gsl/gsl_sf_gamma.h>
#include "util.h"

static int initialized = 0;

static double *x, *dEdx_rad, *dEdx, *csda_range;
static size_t size;

static gsl_interp_accel *acc_dEdx_rad;
static gsl_spline *spline_dEdx_rad;

static gsl_interp_accel *acc_dEdx;
static gsl_spline *spline_dEdx;

static gsl_interp_accel *acc_range;
static gsl_spline *spline_range;

/* Electron critical energy in H2O and D2O.
 *
 * These values come from
 * http://pdgprod.lbl.gov/~deg/AtomicNuclearProperties/HTML/deuterium_oxide_liquid.html for D2O and
 * http://pdg.lbl.gov/2018/AtomicNuclearProperties/HTML/water_liquid.html for H2O.
 */
static const double ELECTRON_CRITICAL_ENERGY_H2O = 78.33;
static const double ELECTRON_CRITICAL_ENERGY_D2O = 78.33;

/* Returns the average number of Cerenkov photons in the range 200-800 nm
 * produced by secondary particles in an electron shower.
 *
 * This comes from fitting the ratio # shower photons/rad loss to the function:
 *
 *     c0 + c1/log(T0*c2 + c3)
 *
 * I don't really have any good theoretical motivation for this. My initial
 * thought was that the number of photons should be roughly proportional to the
 * energy lost due to radiation which is why I chose to fit the ratio. This
 * ratio turned out to vary from approximately 520 at low energies (10 MeV) to
 * a roughly constant value of 420 at higher energies (> 1 GeV).
 *
 * This functional form just happened to fit the ratio as a function of energy
 * pretty well.
 *
 * `T0` is the initial kinetic energy of the electron in MeV and `rad` is the
 * energy lost due to radiation in MeV. */
double electron_get_shower_photons(double T0, double rad)
{
    return rad*(406.745 + 0.271816/log(5.31309e-05*T0+1.00174));
}

void electron_get_position_distribution_parameters(double T0, double *a, double *b)
{
    /* Computes the gamma distribution parameters describing the longitudinal
     * profile of the photons generated from an electromagnetic shower.
     *
     * From the PDG "Passage of Particles" section 32.5:
     *
     * "The mean longitudinal profile of the energy deposition in an
     * electromagnetic cascade is reasonably well described by a gamma
     * distribution."
     *
     * Here we use a slightly different form of the gamma distribution:
     *
     *     f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a)
     *
     * I determined the b parameter by simulating high energy electrons using
     * RAT-PAC and determined that it's roughly equal to the radiation length.
     * To calculate the a parameter we use the formula from the PDG, i.e.
     *
     *     tmax = (a-1)/b = ln(E/E_C) - 0.5
     *
     * Therefore, we calculate a as:
     *
     *     a = tmax*b+1.
     *
     * `T` should be in units of MeV.
     *
     * Example:
     *
     *     double a, b;
     *     electron_get_position_distribution_parameters(1000.0, &a, &b);
     *     double pdf = gamma_pdf(x, a, b);
     *
     * See http://pdg.lbl.gov/2014/reviews/rpp2014-rev-passage-particles-matter.pdf. */
    double tmax;

    *b = RADIATION_LENGTH;
    tmax = log(T0/ELECTRON_CRITICAL_ENERGY_D2O) - 0.5;
    *a = fmax(1.1,tmax*(*b)/RADIATION_LENGTH + 1);
}

double electron_get_angular_distribution_alpha(double T0)
{
    /* To describe the angular distribution of photons in an electromagnetic
     * shower I came up with the heuristic form:
     *
     *     f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
     *
     * I simulated high energy electrons using RAT-PAC in heavy water to fit
     * for the alpha and beta parameters as a function of energy and determined
     * that a reasonably good fit to the data was:
     *
     *     alpha = c0 + c1/log(c2*T0 + c3)
     *
     * where T0 is the initial energy of the electron in MeV and c0, c1, c2,
     * and c3 are constants which I fit out. */
    return fmin(100.0,3.141318e-1 + 2.08198e-01/log(6.33331e-03*T0 + 1.19213e+00));
}

double electron_get_angular_distribution_beta(double T0)
{
    /* To describe the angular distribution of photons in an electromagnetic
     * shower I came up with the heuristic form:
     *
     * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
     *
     * I simulated high energy electrons using RAT-PAC in heavy water to fit
     * for the alpha and beta parameters as a function of energy and determined
     * that a reasonably good fit to the data was:
     *
     * beta = c0 + c1/log(c2*T0 + c3)
     *
     * where T0 is the initial energy of the electron in MeV and c0, c1, c2,
     * and c3 are constants which I fit out. */
    return 1.35372e-01 + 2.22344e-01/log(1.96249e-02*T0 + 1.23912e+00);
}

double electron_get_angular_pdf_no_norm(double cos_theta, double alpha, double beta, double mu)
{
    return exp(-pow(fabs(cos_theta-mu),alpha)/beta);
}

double electron_get_angular_pdf_norm(double alpha, double beta, double mu)
{
    return pow(beta,1/alpha)*gsl_sf_gamma(1/alpha)*(gsl_sf_gamma_inc_P(1/alpha,pow(1-mu,alpha)/beta)+gsl_sf_gamma_inc_P(1/alpha,pow(1+mu,alpha)/beta))/alpha;
}

double electron_get_angular_pdf_delta_ray(double cos_theta, double alpha, double beta, double mu)
{
    /* Returns the probability density that a photon in the wavelength range
     * 200 nm - 800 nm is emitted at an angle cos_theta.
     *
     * The angular distribution is modelled by the pdf:
     *
     *     f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
     *
     * where alpha and beta are constants which depend on the initial energy of
     * the particle.
     *
     * There is no nice closed form solution for the integral of this function,
     * so we just compute it on the fly. To make things faster, we keep track
     * of the last alpha, beta, and mu parameters that were passed in so we
     * avoid recomputing the normalization factor. */
    size_t i;
    static double last_alpha, last_beta, last_mu, norm;
    static int first = 1;
    static double x[10000], y[10000];

    if (first || alpha != last_alpha || beta != last_beta || mu != last_mu) {
        norm = electron_get_angular_pdf_norm(alpha, beta, mu);

        last_alpha = alpha;
        last_beta = beta;
        last_mu = mu;

        for (i = 0; i < LEN(x); i++) {
            x[i] = -1.0 + 2.0*i/(LEN(x)-1);
            y[i] = electron_get_angular_pdf_no_norm(x[i], alpha, beta, mu)/norm;
        }

        first = 0;
    }

    return interp1d(cos_theta,x,y,LEN(x));
}

double electron_get_angular_pdf(double cos_theta, double alpha, double beta, double mu)
{
    /* Returns the probability density that a photon in the wavelength range
     * 200 nm - 800 nm is emitted at an angle cos_theta.
     *
     * The angular distribution is modelled by the pdf:
     *
     *     f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
     *
     * where alpha and beta are constants which depend on the initial energy of
     * the particle.
     *
     * There is no nice closed form solution for the integral of this function,
     * so we just compute it on the fly. To make things faster, we keep track
     * of the last alpha, beta, and mu parameters that were passed in so we
     * avoid recomputing the normalization factor. */
    size_t i;
    static double last_alpha, last_beta, last_mu, norm;
    static int first = 1;
    static double x[10000], y[10000];

    if (first || alpha != last_alpha || beta != last_beta || mu != last_mu) {
        norm = electron_get_angular_pdf_norm(alpha, beta, mu);

        last_alpha = alpha;
        last_beta = beta;
        last_mu = mu;

        for (i = 0; i < LEN(x); i++) {
            x[i] = -1.0 + 2.0*i/(LEN(x)-1);
            y[i] = electron_get_angular_pdf_no_norm(x[i], alpha, beta, mu)/norm;
        }

        first = 0;
    }

    return interp1d(cos_theta,x,y,LEN(x));
}

static int init()
{
    int i, j;
    char line[256];
    char *str;
    double value;
    int n;

    FILE *f = open_file("e_water_liquid.txt", "r");

    if (!f) {
        fprintf(stderr, "failed to open e_water_liquid.txt: %s\n", strerror(errno));
        return -1;
    }

    i = 0;
    n = 0;
    /* For the first pass, we just count how many values there are. */
    while (fgets(line, sizeof(line), f)) {
        size_t len = strlen(line);
        if (len && (line[len-1] != '\n')) {
            fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line);
            goto err;
        }

        i += 1;

        /* Skip the first 8 lines since it's just a header. */
        if (i <= 8) continue;

        if (!len) continue;
        else if (line[0] == '#') continue;

        str = strtok(line," \n");

        while (str) {
            value = strtod(str, NULL);
            str = strtok(NULL," \n");
        }

        n += 1;
    }

    x = malloc(sizeof(double)*n);
    dEdx_rad = malloc(sizeof(double)*n);
    dEdx = malloc(sizeof(double)*n);
    csda_range = malloc(sizeof(double)*n);
    size = n;

    i = 0;
    n = 0;
    /* Now, we actually store the values. */
    rewind(f);
    while (fgets(line, sizeof(line), f)) {
        size_t len = strlen(line);
        if (len && (line[len-1] != '\n')) {
            fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line);
            goto err;
        }

        i += 1;

        /* Skip the first 8 lines since it's just a header. */
        if (i <= 8) continue;

        if (!len) continue;
        else if (line[0] == '#') continue;

        str = strtok(line," \n");

        j = 0;
        while (str) {
            value = strtod(str, NULL);
            switch (j) {
            case 0:
                x[n] = value;
                break;
            case 2:
                dEdx_rad[n] = value;
                break;
            case 3:
                dEdx[n] = value;
                break;
            case 4:
                csda_range[n] = value;
                break;
            }
            j += 1;
            str = strtok(NULL," \n");
        }

        n += 1;
    }

    fclose(f);

    acc_dEdx_rad = gsl_interp_accel_alloc();
    spline_dEdx_rad = gsl_spline_alloc(gsl_interp_linear, size);
    gsl_spline_init(spline_dEdx_rad, x, dEdx_rad, size);

    acc_dEdx = gsl_interp_accel_alloc();
    spline_dEdx = gsl_spline_alloc(gsl_interp_linear, size);
    gsl_spline_init(spline_dEdx, x, dEdx, size);

    acc_range = gsl_interp_accel_alloc();
    spline_range = gsl_spline_alloc(gsl_interp_linear, size);
    gsl_spline_init(spline_range, x, csda_range, size);

    initialized = 1;

    return 0;

err:
    fclose(f);

    return -1;
}

/* Returns the maximum kinetic energy for an electron in the range tables.
 *
 * If you call electron_get_range() or electron_get_dEdx() with a kinetic
 * energy higher you will get a GSL interpolation error. */
double electron_get_max_energy(void)
{
    if (!initialized) {
        if (init()) {
            exit(1);
        }
    }

    return x[size-1];
}

double electron_get_range(double T, double rho)
{
    /* Returns the approximate range a electron with kinetic energy `T` will travel
     * in water before losing all of its energy. This range is interpolated
     * based on data from the PDG which uses the continuous slowing down
     * approximation.
     *
     * `T` should be in MeV, and `rho` should be in g/cm^3.
     *
     * Return value is in cm.
     *
     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
    if (!initialized) {
        if (init()) {
            exit(1);
        }
    }

    return gsl_spline_eval(spline_range, T, acc_range)/rho;
}

double electron_get_dEdx_rad(double T, double rho)
{
    /* Returns the approximate radiative dE/dx for a electron in water with
     * kinetic energy `T`.
     *
     * `T` should be in MeV and `rho` in g/cm^3.
     *
     * Return value is in MeV/cm.
     *
     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
    if (!initialized) {
        if (init()) {
            exit(1);
        }
    }

    if (T < spline_dEdx_rad->x[0]) return spline_dEdx_rad->y[0];

    return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho;
}

double electron_get_dEdx(double T, double rho)
{
    /* Returns the approximate dE/dx for a electron in water with kinetic energy
     * `T`.
     *
     * `T` should be in MeV and `rho` in g/cm^3.
     *
     * Return value is in MeV/cm.
     *
     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
    if (!initialized) {
        if (init()) {
            exit(1);
        }
    }

    if (T < spline_dEdx->x[0]) return spline_dEdx->y[0];

    return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
}