#include "Record_Info.h" #include #include "zdab_utils.h" #include "pack2b.h" #include /* for size_t */ #include /* for fprintf() */ int isOrphan(aPmtEventRecord *pmtRecord) { /* Returns non-zero if the specified event is an orphan. */ int i; uint32_t *mtc_data = (uint32_t *) &pmtRecord->TriggerCardData; for (i = 0; i < 6; ++i) { if (*mtc_data != 0) return 0; ++mtc_data; } return 1; } // PH 04/23/98 // Swap 4-byte integer/floats between native and external format void swap_int32(int32_t *val_pt, int count) { int32_t *last = val_pt + count; while (val_pt < last) { *val_pt = ((*val_pt << 24) & 0xff000000) | ((*val_pt << 8) & 0x00ff0000) | ((*val_pt >> 8) & 0x0000ff00) | ((*val_pt >> 24) & 0x000000ff); ++val_pt; } return; } // Swap 2-byte integers between native and external format void swap_int16(int16_t *val_pt, int count) { char tmp; int i; for (i=0; icls); unpack((uint8_t *) (data+1), "l",&b->inc); unpack((uint8_t *) (data+2), "f",&b->x); unpack((uint8_t *) (data+3), "f",&b->y); unpack((uint8_t *) (data+4), "f",&b->z); unpack((uint8_t *) (data+5), "F",&b->tim); unpack((uint8_t *) (data+7), "l",&b->rgn); unpack((uint8_t *) (data+8), "l",&b->idm); unpack((uint8_t *) (data+9), "l",&b->rg2); unpack((uint8_t *) (data+10), "l",&b->im2); unpack((uint8_t *) (data+12), "f",&b->bnx); unpack((uint8_t *) (data+13), "f",&b->bny); unpack((uint8_t *) (data+14), "f",&b->bnz); unpack((uint8_t *) (data+15), "l",&b->cer); } void unpack_mctk(uint32_t *data, MCTKBank *b) { unpack((uint8_t *) data, "l",&b->idp); unpack((uint8_t *) (data+1), "f",&b->drx); unpack((uint8_t *) (data+2), "f",&b->dry); unpack((uint8_t *) (data+3), "f",&b->drz); unpack((uint8_t *) (data+4), "f",&b->ene); unpack((uint8_t *) (data+5), "l",&b->rgn); unpack((uint8_t *) (data+6), "l",&b->idm); unpack((uint8_t *) (data+7), "f",&b->plx); unpack((uint8_t *) (data+8), "f",&b->ply); unpack((uint8_t *) (data+9), "f",&b->plz); unpack((uint8_t *) (data+10), "f",&b->stp); unpack((uint8_t *) (data+11), "f",&b->near); } void unpack_ev(uint32_t *data, EVBank *b) { unpack((uint8_t *) data, "l",&b->run); unpack((uint8_t *) (data+1), "l",&b->evn); unpack((uint8_t *) (data+2), "l",&b->dtp); unpack((uint8_t *) (data+3), "l",&b->jdy); unpack((uint8_t *) (data+4), "l",&b->ut1); unpack((uint8_t *) (data+5), "l",&b->ut2); unpack((uint8_t *) (data+6), "l",&b->dte); unpack((uint8_t *) (data+7), "l",&b->hmsc); unpack((uint8_t *) (data+8), "l",&b->gtr1); unpack((uint8_t *) (data+9), "l",&b->gtr2); unpack((uint8_t *) (data+10),"l",&b->npm); unpack((uint8_t *) (data+11),"l",&b->nph); unpack((uint8_t *) (data+12),"l",&b->sub_run); unpack((uint8_t *) (data+13),"l",&b->mc_pck); unpack((uint8_t *) (data+14),"l",&b->rec); unpack((uint8_t *) (data+15),"l",&b->vpck); unpack((uint8_t *) (data+16),"l",&b->gtr_id); unpack((uint8_t *) (data+17),"l",&b->trg_type); unpack((uint8_t *) (data+18),"l",&b->peak); unpack((uint8_t *) (data+19),"l",&b->diff); unpack((uint8_t *) (data+20),"l",&b->integral); unpack((uint8_t *) (data+21),"l",&b->err); unpack((uint8_t *) (data+22),"l",&b->data_set); unpack((uint8_t *) (data+22),"lll",&b->spare1[0], &b->spare1[1], &b->spare1[2]); unpack((uint8_t *) (data+26),"l",&b->ncd_status); unpack((uint8_t *) (data+27),"l",&b->num_muxg); unpack((uint8_t *) (data+29),"l",&b->num_mux); unpack((uint8_t *) (data+29),"l",&b->num_scope); unpack((uint8_t *) (data+30),"lllll",&b->spare2[0], &b->spare2[1], &b->spare2[2], &b->spare2[3], &b->spare2[4]); unpack((uint8_t *) (data+35),"l",&b->ncd_clk_up); unpack((uint8_t *) (data+36),"l",&b->ncd_clk_lw); unpack((uint8_t *) (data+37),"l",&b->ncd_reg); unpack((uint8_t *) (data+38),"l",&b->ncd_gtid); unpack((uint8_t *) (data+39),"l",&b->ncd_sync); unpack((uint8_t *) (data+40),"l",&b->spare3[0], &b->spare3[1], &b->spare3[2], &b->spare3[3], &b->spare3[4], &b->spare3[5], &b->spare3[6], &b->spare3[7], &b->spare3[8], &b->spare3[9]); } void unpack_pmt(uint32_t *data, PMTBank *b) { unpack((uint8_t *) data,"l",&b->pn); unpack((uint8_t *) (data+1),"l",&b->pf); unpack((uint8_t *) (data+2),"f",&b->pt); unpack((uint8_t *) (data+3),"f",&b->phl); unpack((uint8_t *) (data+4),"f",&b->phs); unpack((uint8_t *) (data+5),"f",&b->plx); unpack((uint8_t *) (data+6),"f",&b->pt0); unpack((uint8_t *) (data+7),"l",&b->pif); unpack((uint8_t *) (data+8),"f",&b->pit); unpack((uint8_t *) (data+9),"f",&b->pihl); unpack((uint8_t *) (data+10),"f",&b->pihs); unpack((uint8_t *) (data+11),"f",&b->pilx); unpack((uint8_t *) (data+12),"f",&b->pit0); unpack((uint8_t *) (data+13),"l",&b->cell); unpack((uint8_t *) (data+14),"l",&b->pin); unpack((uint8_t *) (data+15),"l",&b->tslh); unpack((uint8_t *) (data+16),"l",&b->hca); unpack((uint8_t *) (data+17),"l",&b->eca_val); unpack((uint8_t *) (data+18),"l",&b->pca_val); unpack((uint8_t *) (data+19),"l",&b->anxx); unpack((uint8_t *) (data+20),"l",&b->ept); unpack((uint8_t *) (data+21),"l",&b->ehl); unpack((uint8_t *) (data+22),"l",&b->ehs); unpack((uint8_t *) (data+23),"l",&b->elx); unpack((uint8_t *) (data+24),"l",&b->pt1); unpack((uint8_t *) (data+25),"l",&b->ptm); unpack((uint8_t *) (data+26),"l",&b->ptms); unpack((uint8_t *) (data+27),"l",&b->qm); unpack((uint8_t *) (data+28),"l",&b->qms); unpack((uint8_t *) (data+29),"l",&b->qrc); } int swap_PmtRecord(aPmtEventRecord *aPmtRecord, size_t size) { /* Swap a Pmt Event Record. This function swaps both the Pmt event record * and the PMT hits and sub fields. Returns -1 if the PMT record has too * many hits. */ SWAP_INT32(aPmtRecord, sizeof(aPmtEventRecord)/sizeof(uint32_t)); int npmt = aPmtRecord->NPmtHit; if (npmt > MAX_NHIT) { fprintf(stderr, "Read error: Bad ZDAB -- %d pmt hit!", npmt); return -1; } else { if (size < sizeof(aPmtEventRecord) + 3*npmt*4) { fprintf(stderr, "swap_PmtRecord: size of record is %zu bytes, but there are %i PMT hits", size, npmt); return -1; } // swap the hit data SWAP_INT32(aPmtRecord + 1, 3*npmt); // swap the sub-fields uint32_t *sub_header = &aPmtRecord->CalPckType; while (*sub_header & SUB_NOT_LAST) { if (size < (sub_header - (uint32_t *) aPmtRecord)*4 + (*sub_header & SUB_LENGTH_MASK)*4 + 4) { fprintf(stderr, "swap_PmtRecord: size of record is %zu bytes, " "but sub-field requires %lu bytes", size, (sub_header - (uint32_t *) aPmtRecord)*4 + (*sub_header & SUB_LENGTH_MASK)*4 + 4); return -1; } sub_header += (*sub_header & SUB_LENGTH_MASK); SWAP_INT32(sub_header, 1); // swap the sub-field header // get number of data words (-1 because we don't want to include header size) uint32_t data_words = (*sub_header & SUB_LENGTH_MASK) - 1; if (size < (sub_header - (uint32_t *) aPmtRecord)*4 + (*sub_header & SUB_LENGTH_MASK)*4) { fprintf(stderr, "swap_PmtRecord: size of record is %zu bytes, " "but sub-field requires %lu bytes", size, (sub_header - (uint32_t *) aPmtRecord)*4 + (*sub_header & SUB_LENGTH_MASK)*4); return -1; } SWAP_INT32(sub_header+1, data_words); } } return 0; } void swap_TrigRecord(struct TriggerInfo *aTrigRecord) { /* Byte swap a Trigger Record. */ SWAP_INT32(aTrigRecord, sizeof(struct TriggerInfo)/sizeof(uint32_t)); } void swap_RunRecord(struct RunRecord *aRunRecord) { /* Byte swap a Run Record. */ SWAP_INT32(aRunRecord, sizeof(struct RunRecord)/sizeof(uint32_t)); } >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 
#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);
}
generated by cgit (git 2.34.1) at 2025-05-08 15:23:13 +0000