src/pmt_response.c


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83Building a Geometry
===================

.. module:: chroma.geometry

The Mesh Class
--------------

.. autoclass:: Mesh
   :members:

The Solid Class
---------------

.. autoclass:: Solid
   :members:

The Material Class
------------------

.. autoclass:: Material
   :members:

The Geometry Class
------------------

.. autoclass:: Geometry
   :members:

Importing from STL
------------------

.. autofunction:: chroma.mesh_from_stl

Mesh Modelling Tools
--------------------

.. module:: chroma.make

.. autofunction:: linear_extrude

.. image:: images/hexagonal_prism.png
   :height: 100px

.. autofunction:: rotate_extrude

.. image:: images/bipyramid.png
   :height: 100px

.. autofunction:: box

.. image:: images/box.png
   :height: 100px

.. autofunction:: cube

.. image:: images/cube.png
   :height: 100px

.. autofunction:: cylinder

.. image:: images/cylinder.png
   :height: 100px

.. autofunction:: segmented_cylinder

.. autofunction:: sphere

.. image:: images/sphere.png
   :height: 100px

.. autofunction:: torus

.. image:: images/torus.png
   :height: 100px
#include "pmt_response.h"
#include "dict.h"
#include <stdint.h> /* for uint32_t */
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_integration.h>
#include <stdio.h> /* for sprintf() */
#include "misc.h"
#include "quantum_efficiency.h"
#include "db.h"
#include <gsl/gsl_errno.h> /* for gsl_strerror() */

static int initialized = 0;

char pmtr_err[256];

static gsl_spline2d *spline_resp;
static gsl_spline2d *spline_reflec;
static gsl_interp_accel *xacc;
static gsl_interp_accel *yacc;

static double thetas[NUM_ANGLES] = {
    0.0,
    1.0,
    2.0,
    3.0,
    4.0,
    5.0,
    6.0,
    7.0,
    8.0,
    9.0,
    10.0,
    11.0,
    12.0,
    13.0,
    14.0,
    15.0,
    16.0,
    17.0,
    18.0,
    19.0,
    20.0,
    21.0,
    22.0,
    23.0,
    24.0,
    25.0,
    26.0,
    27.0,
    28.0,
    29.0,
    30.0,
    31.0,
    32.0,
    33.0,
    34.0,
    35.0,
    36.0,
    37.0,
    38.0,
    39.0,
    40.0,
    41.0,
    42.0,
    43.0,
    44.0,
    45.0,
    46.0,
    47.0,
    48.0,
    49.0,
    50.0,
    51.0,
    52.0,
    53.0,
    54.0,
    55.0,
    56.0,
    57.0,
    58.0,
    59.0,
    60.0,
    61.0,
    62.0,
    63.0,
    64.0,
    65.0,
    66.0,
    67.0,
    68.0,
    69.0,
    70.0,
    71.0,
    72.0,
    73.0,
    74.0,
    75.0,
    76.0,
    77.0,
    78.0,
    79.0,
    80.0,
    81.0,
    82.0,
    83.0,
    84.0,
    85.0,
    86.0,
    87.0,
    88.0,
    89.0,
};

static double wavelengths[NUM_WAVELENGTHS] = {
    220.0,
    230.0,
    240.0,
    250.0,
    260.0,
    270.0,
    280.0,
    290.0,
    300.0,
    310.0,
    320.0,
    330.0,
    340.0,
    350.0,
    360.0,
    370.0,
    380.0,
    390.0,
    400.0,
    410.0,
    420.0,
    430.0,
    440.0,
    450.0,
    460.0,
    470.0,
    480.0,
    490.0,
    500.0,
    510.0,
    520.0,
    530.0,
    540.0,
    550.0,
    560.0,
    570.0,
    580.0,
    590.0,
    600.0,
    610.0,
    620.0,
    630.0,
    640.0,
    650.0,
    660.0,
    670.0,
    680.0,
    690.0,
    700.0,
    710.0,
};

/* PMT response as a function of angle and wavelength. */
static double resp[NUM_ANGLES*NUM_WAVELENGTHS];
/* PMT response as a function of angle weighted by the Cerenkov spectrum and
 * the quantum efficiency. */
static double weighted_resp[NUM_ANGLES];

/* PMT reflectivity as a function of angle and wavelength. */
static double reflec[NUM_ANGLES*NUM_WAVELENGTHS];
/* PMT reflectivity as a function of angle weighted by the Cerenkov spectrum
 * and the quantum efficiency. */
static double weighted_reflec[NUM_ANGLES];

double get_weighted_pmt_reflectivity(double theta)
{
    /* Returns the probability that a photon is reflected as a function of the
     * photon angle with respect to the PMT normal. `theta` should be in
     * radians. */
    double deg;

    if (!initialized) {
        fprintf(stderr, "pmt response hasn't been initialized!\n");
        exit(1);
    }

    /* Convert to degrees since the PMTR table is in degrees. */
    deg = theta*180.0/M_PI;

    deg = fmod(deg,180.0);

    if (deg < 0.0)
        deg += 180.0;

    return interp1d(deg, thetas, weighted_reflec, NUM_ANGLES);
}

double get_weighted_pmt_response(double theta)
{
    /* Returns the probability that a photon is detected as a function of the
     * photon angle with respect to the PMT normal. `theta` should be in
     * radians.
     *
     * Note: This does *not* include the PMT quantum efficiency. */
    double deg;

    if (!initialized) {
        fprintf(stderr, "pmt response hasn't been initialized!\n");
        exit(1);
    }

    /* Convert to degrees since the PMTR table is in degrees. */
    deg = theta*180.0/M_PI;

    deg = fmod(deg,180.0);

    if (deg < 0.0)
        deg += 180.0;

    return interp1d(deg, thetas, weighted_resp, NUM_ANGLES);
}

double get_pmt_reflectivity(double wavelength, double theta)
{
    /* Returns the probability that a photon is reflected as a function of the
     * photon angle with respect to the PMT normal (in radians) and the
     * wavelength (in nm). */
    double deg;

    if (!initialized) {
        fprintf(stderr, "pmt response hasn't been initialized!\n");
        exit(1);
    }

    /* Convert to degrees since the PMTR table is in degrees. */
    deg = theta*180.0/M_PI;

    deg = fmod(deg,180.0);

    if (deg < 0.0)
        deg += 180.0;

    if (deg > thetas[NUM_ANGLES-1])
        deg = thetas[NUM_ANGLES-1];

    if (wavelength < wavelengths[0])
        wavelength = wavelengths[0];

    if (wavelength > wavelengths[NUM_WAVELENGTHS-1])
        wavelength = wavelengths[NUM_WAVELENGTHS-1];

    return gsl_spline2d_eval(spline_reflec, deg, wavelength, xacc, yacc);
}

double get_pmt_response(double wavelength, double theta)
{
    /* Returns the probability that a photon is detected as a function of the
     * photon angle with respect to the PMT normal (in radians) and the
     * wavelength (in nm).
     *
     * Note: This does *not* include the PMT quantum efficiency. */
    double deg, qe;

    if (!initialized) {
        fprintf(stderr, "pmt response hasn't been initialized!\n");
        exit(1);
    }

    /* Convert to degrees since the PMTR table is in degrees. */
    deg = theta*180.0/M_PI;

    deg = fmod(deg,180.0);

    if (deg < 0.0)
        deg += 180.0;

    if (deg > thetas[NUM_ANGLES-1])
        deg = thetas[NUM_ANGLES-1];

    if (wavelength < wavelengths[0])
        wavelength = wavelengths[0];

    if (wavelength > wavelengths[NUM_WAVELENGTHS-1])
        wavelength = wavelengths[NUM_WAVELENGTHS-1];

    /* The PMTR bank in SNOMAN included the effect of the PMT quantum
     * efficiency since it was used for the "grey disk" model.  Therefore,
     * since we already account for the quantum efficiency it is necessary to
     * divide the response by the quantum efficiency. */

    qe = get_quantum_efficiency(wavelength);

    /* If the quantum efficiency is zero, we have no way of knowing what the
     * response should be since it was multiplied by zero. Since for these
     * wavelengths the photon won't be detected, it doesn't really matter. */
    if (qe == 0.0) return 0.0;

    return gsl_spline2d_eval(spline_resp, deg, wavelength, xacc, yacc)/qe;
}

static double gsl_pmt_reflectivity(double wavelength, void *params)
{
    /* Returns the probability that a photon is reflected as a function of
     * wavelength for a specific angle weighted by the quantum efficiency and
     * the Cerenkov spectrum. */
    double qe, theta;

    theta = ((double *) params)[0];

    qe = get_quantum_efficiency(wavelength);

    return qe*get_pmt_reflectivity(wavelength,theta)/pow(wavelength,2);
}

static double gsl_pmt_response(double wavelength, void *params)
{
    /* Returns the probability that a photon is detected as a function of
     * wavelength for a specific angle weighted by the quantum efficiency and
     * the Cerenkov spectrum. */
    double qe, theta;

    theta = ((double *) params)[0];

    qe = get_quantum_efficiency(wavelength);

    return qe*get_pmt_response(wavelength,theta)/pow(wavelength,2);
}

static double gsl_cerenkov(double wavelength, void *params)
{
    /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */
    double qe;

    qe = get_quantum_efficiency(wavelength);

    return qe/pow(wavelength,2);
}

int pmt_response_init(dict *db)
{
    int i, j;
    float *pmtr;
    double norm;

    double result, error;
    size_t nevals;
    int status;
    gsl_integration_cquad_workspace *w;
    gsl_function F;
    double params[1];

    spline_resp = gsl_spline2d_alloc(gsl_interp2d_bilinear, NUM_ANGLES, NUM_WAVELENGTHS);
    spline_reflec = gsl_spline2d_alloc(gsl_interp2d_bilinear, NUM_ANGLES, NUM_WAVELENGTHS);
    xacc = gsl_interp_accel_alloc();
    yacc = gsl_interp_accel_alloc();

    pmtr = (float *) get_bank(db, "PMTR", 1);

    if (!pmtr) {
        sprintf(pmtr_err, "failed to load PMTR bank\n");
        return -1;
    }

    for (i = 0; i < NUM_ANGLES; i++) {
        for (j = 0; j < NUM_WAVELENGTHS; j++) {
            gsl_spline2d_set(spline_resp, resp, i, j, pmtr[30+KPMTR_RESP+i+j*NUM_ANGLES]);
            gsl_spline2d_set(spline_reflec, reflec, i, j, pmtr[30+KPMTR_REFLEC+i+j*NUM_ANGLES]);
        }
    }

    gsl_spline2d_init(spline_resp, thetas, wavelengths, resp, NUM_ANGLES, NUM_WAVELENGTHS);
    gsl_spline2d_init(spline_reflec, thetas, wavelengths, reflec, NUM_ANGLES, NUM_WAVELENGTHS);

    initialized = 1;

    w = gsl_integration_cquad_workspace_alloc(100);

    F.function = &gsl_cerenkov;
    F.params = params;

    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &norm, &error, &nevals);

    if (status) {
        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
        exit(1);
    }

    F.function = &gsl_pmt_response;

    for (i = 0; i < NUM_ANGLES; i++) {
        params[0] = thetas[i]*M_PI/180.0;

        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);

        weighted_resp[i] = result/norm;

        if (status) {
            fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
            exit(1);
        }
    }

    F.function = &gsl_pmt_reflectivity;

    for (i = 0; i < NUM_ANGLES; i++) {
        params[0] = thetas[i]*M_PI/180.0;

        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);

        weighted_reflec[i] = result/norm;

        if (status) {
            fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
            exit(1);
        }
    }

    gsl_integration_cquad_workspace_free(w);

    return 0;
}

void pmt_response_free(void)
{
    gsl_spline2d_free(spline_resp);
    gsl_spline2d_free(spline_reflec);
    gsl_interp_accel_free(xacc);
    gsl_interp_accel_free(yacc);
}