aboutsummaryrefslogtreecommitdiff
path: root/src/random.c
blob: e2e56418e6cbee831ac25d7f0a424b1b4bd4b249 (plain)
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
#include "random.h"
#include <math.h>

double randn(void)
{
    /* Generates a random number from a normal distribution using the
     * Box-Muller transform. */
    double u1, u2;

    u1 = genrand_real1();
    u2 = genrand_real1();

    return sqrt(-2*log(u1))*cos(2*M_PI*u2);
}

void rand_sphere(double *dir)
{
    /* Generates a random point on the unit sphere. */
    double u, v, theta, phi;

    u = genrand_real1();
    v = genrand_real1();

    phi = 2*M_PI*u;
    theta = acos(2*v-1);

    dir[0] = sin(theta)*cos(phi);
    dir[1] = sin(theta)*sin(phi);
    dir[2] = cos(theta);
}
es { font-weight: bold; font-style: italic } /* Generic.EmphStrong */ .highlight .gr { color: #aa0000 } /* Generic.Error */ .highlight .gh { color: #333333 } /* Generic.Heading */ .highlight .gi { color: #000000; background-color: #ddffdd } /* Generic.Inserted */ .highlight .go { color: #888888 } /* Generic.Output */ .highlight .gp { color: #555555 } /* Generic.Prompt */ .highlight .gs { font-weight: bold } /* Generic.Strong */ .highlight .gu { color: #666666 } /* Generic.Subheading */ .highlight .gt { color: #aa0000 } /* Generic.Traceback */ .highlight .kc { color: #008800; font-weight: bold } /* Keyword.Constant */ .highlight .kd { color: #008800; font-weight: bold } /* Keyword.Declaration */ .highlight .kn { color: #008800; font-weight: bold } /* Keyword.Namespace */ .highlight .kp { color: #008800 } /* Keyword.Pseudo */ .highlight .kr { color: #008800; font-weight: bold } /* Keyword.Reserved */ .highlight .kt { color: #888888; font-weight: bold } /* Keyword.Type */ .highlight .m { color: #0000DD; font-weight: bold } /* Literal.Number */ .highlight .s { color: #dd2200; background-color: #fff0f0 } /* Literal.String */ .highlight .na { color: #336699 } /* Name.Attribute */ .highlight .nb { color: #003388 } /* Name.Builtin */ .highlight .nc { color: #bb0066; font-weight: bold } /* Name.Class */ .highlight .no { color: #003366; font-weight: bold } /* Name.Constant */ .highlight .nd { color: #555555 } /* Name.Decorator */ .highlight .ne { color: #bb0066; font-weight: bold } /* Name.Exception */ .highlight .nf { color: #0066bb; font-weight: bold } /* Name.Function */ .highlight .nl { color: #336699; font-style: italic } /* Name.Label */ .highlight .nn { color: #bb0066; font-weight: bold } /* Name.Namespace */ .highlight .py { color: #336699; font-weight: bold } /* Name.Property */ .highlight .nt { color: #bb0066; font-weight: bold } /* Name.Tag */ .highlight .nv { color: #336699 } /* Name.Variable */ .highlight .ow { color: #008800 } /* Operator.Word */ .highlight .w { color: #bbbbbb } /* Text.Whitespace */ .highlight .mb { color: #0000DD; font-weight: bold } /* Literal.Number.Bin */ .highlight .mf { color: #0000DD; font-weight: bold } /* Literal.Number.Float */ .highlight .mh { color: #0000DD; font-weight: bold } /* Literal.Number.Hex */ .highlight .mi { color: #0000DD; font-weight: bold } /* Literal.Number.Integer */ .highlight .mo { color: #0000DD; font-weight: bold } /* Literal.Number.Oct */ .highlight .sa { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Affix */ .highlight .sb { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Backtick */ .highlight .sc { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Char */ .highlight .dl { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Delimiter */ .highlight .sd { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Doc */ .highlight .s2 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Double */ .highlight .se { color: #0044dd; background-color: #fff0f0 } /* Literal.String.Escape */ .highlight .sh { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Heredoc */ .highlight .si { color: #3333bb; background-color: #fff0f0 } /* Literal.String.Interpol */ .highlight .sx { color: #22bb22; background-color: #f0fff0 } /* Literal.String.Other */ .highlight .sr { color: #008800; background-color: #fff0ff } /* Literal.String.Regex */ .highlight .s1 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Single */ .highlight .ss { color: #aa6600; background-color: #fff0f0 } /* Literal.String.Symbol */ .highlight .bp { color: #003388 } /* Name.Builtin.Pseudo */ .highlight .fm { color: #0066bb; font-weight: bold } /* Name.Function.Magic */ .highlight .vc { color: #336699 } /* Name.Variable.Class */ .highlight .vg { color: #dd7700 } /* Name.Variable.Global */ .highlight .vi { color: #3333bb } /* Name.Variable.Instance */ .highlight .vm { color: #336699 } /* Name.Variable.Magic */ .highlight .il { color: #0000DD; font-weight: bold } /* Literal.Number.Integer.Long */
#ifndef ZDAB_UTILS
#define ZDAB_UTILS

#include "Record_Info.h"
#include <stdint.h>
#include <stdlib.h> /* for size_t */

// the builder won't put out events with NHIT > 10000
// (note that these are possible due to hardware problems)
// but XSNOED can write an event with up to 10240 channels
#define MAX_NHIT 10240

#ifdef SWAP_BYTES
#define SWAP_INT32(a,b) swap_int32((int32_t *)(a),(b))
#define SWAP_INT16(a,b) swap_int16((int16_t *)(a),(b))
#else
#define SWAP_INT32(a,b)
#define SWAP_INT16(a,b)
#endif

/* Status bits for the MC Vertex bank.
 *
 * See http://www.lip.pt/~barros/SNO/doc/snoman/companion_frames.html. */
#define KMCVX_SRC 0x00000001  /* Source */
#define KMCVX_BOU 0x00000002  /* Boundary */
#define KMCVX_INT 0x00000004  /* Interaction */
#define KMCVX_SNK 0x00000008  /* Sink */
#define KMCVX_PSC 0x00000010  /* Pre-source */
#define KMCVX_CRE 0x00000100  /* Creation i.e. vertex creates new particle (ignoring Cerenkov photon). */
#define KMCVX_IPM 0x00000200  /* Indirect MCPM. The MCPM hit comes indirectly from the vertex. */

/* Status bits for the MC Track bank.
 *
 * Bits 4-17 are used to store the "history" of Cerenkov photons. The bit
 * structure is the same as in the MCPM bank status word. More details are
 * available on the MCPM bank page.
 *
 * See http://www.lip.pt/~barros/SNO/doc/snoman/companion_frames.html. */
#define KMCTK_IVX        0x00000002  /* Indirect vertex. Track does not come directly from its
                                        supporting vertex, see Notes below. */
#define KMCTK_SPARE1     0x00000004  /* Spare (Not part of Cerenkov history) */
#define KMCTK_SPARE2     0x00000008  /* Spare (Not part of Cerenkov history) */
#define KMCTK_INDIRECT   0x00000010  /* Indirect photon. */
#define KMCTK_RSCAT      0x00000020  /* Photon was Raleigh scattered. */
#define KMCTK_PSUP_SREF  0x00000040  /* Photon had a spectral reflection from the PSUP. */
#define KMCTK_PSUP_DREF  0x00000080  /* Photon had a diffuse reflection from the PSUP. */
#define KMCTK_AV_SREF    0x00000100  /* Photon had a spectral reflection from the AV. */
#define KMCTK_AV_DREF    0x00000200  /* Photon had a diffuse reflection from the AV. */
#define KMCTK_NCD_SREF   0x00000400  /* Photon had a spectral reflection from the NCD OVL. */
#define KMCTK_NCD_DREF   0x00000800  /* Photon had a diffuse reflection from the NCD OVL. */
#define KMCTK_PMT_REF    0x00001000  /* Photon entered and then escaped the PMT code (vxpmt). */
#define KMCTK_IN_PMT     0x00002000  /* Photon was produced inside a PMT. */
#define KMCTK_THROUGH_AV 0x00004000  /* Photon passed through the AV at least once. */

/* Bit masks for the PMT flags field in the PMT bank. */
#define KPF_DIS      0x00000001
#define KPF_NOI      0x00000004
#define KPF_PRP      0x00000008
#define KPF_AFP      0x00000010
#define KPF_FECD     0x00000020
#define KPF_FT_MASK1 0x00000040
#define KPF_XT       0x04000000
#define KPF_NO_CAL   0x08000000
#define KPF_BAD_CAL  0x10000000

typedef struct MCVXBank {
    /* Class:  = 1 Source, = 2 Boundary, = 3 Interaction, = 4 Sink
	       = 5 Pre-source. */
    uint32_t cls;
    /* Interaction code. */
    uint32_t inc;
    /* Position X. */
    float x;
    /* Position Y. */
    float y;
    /* Position Z. */
    float z;
    /* Time in ns since MC Generation time (in MC bank). */
    double tim;
    /* First region code. */
    uint32_t rgn;
    /* First physical media code. */
    uint32_t idm;
    /* Second region code. */
    uint32_t rg2;
    /* Second physical media code. */
    uint32_t im2;
    /* Boundary normal X (only defined if boundary status bit set). */
    float bnx;
    /* Boundary normal Y (only defined if boundary status bit set). */
    float bny;
    /* Boundary normal Z (only defined if boundary status bit set). */
    float bnz;
    /* Number of Cerenkov photons. Only non-zero for CBV vertex. */
    uint32_t cer;
} MCVXBank;

typedef struct MCTKBank {
    /* Particle id code. */
    uint32_t idp;
    /* Direction cosine X. */
    float drx;
    /* Direction cosine Y. */
    float dry;
    /* Direction cosine Z. */
    float drz;
    /* Total energy (except for neutrons where it is kinetic). */
    float ene;
    /* Region code. */
    uint32_t rgn;
    /* Physical media code. */
    uint32_t idm;
    /* Polarization X. */
    float plx;
    /* Polarization Y. */
    float ply;
    /* Polarization Z. */
    float plz;
    /* Track step size. */
    float stp;
    /* Minimum distance to nearest boundary or 0. */
    float near;
} MCTKBank;

typedef struct EVBank {
    /* Run number. */
    uint32_t run;
    /* Event number. */
    uint32_t evn;
    /* Data type. */
    uint32_t dtp;
    /* Julian date. */
    uint32_t jdy;
    /* Universal time (seconds). */
    uint32_t ut1;
    /* Universal time (nanoseconds). */
    uint32_t ut2;
    /* Date (format: yyyymmdd) */
    uint32_t dte;
    /* Time (format: hhmmsscc - cc is centisec). */
    uint32_t hmsc;
    /* Global trigger time in nsec (first word). */
    uint32_t gtr1;
    /* Second word. */
    uint32_t gtr2;
    /* Number of PMTs that fired. */
    uint32_t npm;
    /* Integrated charge in the event. */
    uint32_t nph;
    /* Sub-run number, or -1 if unknown. Only defined for run number >= 10614. */
    uint32_t sub_run;
    /* Packed word. */
    uint32_t mc_pck;
    /* ZDAB input/output record type, e.g. PMT, NCD, CMA, RUN, etc. */
    uint32_t rec;
    /* ZDAB_PMT format number (x10). */
    uint32_t vpck;
    /* Global trigger ID. */
    uint32_t gtr_id;
    /* Trigger word. */
    uint32_t trg_type;
    /* Digitized peak of analog sum. */
    uint32_t peak;
    /* Analog sum derivative. */
    uint32_t diff;
    /* Digitized integral of analog sum. */
    uint32_t integral;
    /* Trigger error bits plus spares. */
    uint32_t err;
    /* Data splitter blindness word. */
    uint32_t data_set;
    uint32_t spare1[3];
    uint32_t ncd_status;
    /* Number of multiplexer global records. */
    uint32_t num_muxg;
    /* Number of multiplexer records. */
    uint32_t num_mux;
    /* Number of scope records. */
    uint32_t num_scope;
    uint32_t spare2[5];
    /* Upper 24 bits of flock */
    uint32_t ncd_clk_up;
    /* Lower clock register (32 bits). */
    uint32_t ncd_clk_lw;
    /* Latch register id. */
    uint32_t ncd_reg;
    /* GTID. */
    uint32_t ncd_gtid;
    /* Sync clear error. */
    uint32_t ncd_sync;
    uint32_t spare3[10];
} EVBank;

/* PMT Bank struct. This is a struct which mimics the PMT data bank in SNOMAN.
 * The fields and comments were taken from the SNOMAN companion. */
typedef struct PMTBank {
    /* Tube number. For FECD, which has no associated tube number, this word is
     * a copy of PIN. */
    uint32_t pn;
    /* Set of 1-bit flags. All undefined bits are set to zero. The following
     * bits are defined:
     *
     * 0  KPF_DIS       Discard hit for all processing including fitting.
     * 2  KPF_NOI       Noise.
     * 3  KPF_PRP       Pre-pulse.
     * 4  KPF_AFP       After-pulse.
     * 5  KPF_FECD      Front-End Card Data i.e. non-PMT channel.
     *                  Should only be set on the FECD bank.
     * 6  KPF_FT_MASK+i Discarded by fitter method i (i = 1,2,..).
     *                  For description of fitter methods see the FTx bank.
     * 26 KPF_XT        Crosstalk.
     * 27 KPF_NO_CAL    No calibration.
     * 28 KPF_BAD_CAL   Bad calibration. */
    uint32_t pf;
    /* Time in nano-secs relative to event T0 - see EV. */
    float pt;
    /* Integrated charge. */
    float phl;
    /* Short-time integrated charge. */
    float phs;
    /* Low-gain integrated charge. */
    float plx;
    /* Time with the PMT jitter removed. */
    float pt0;
    /* Set of n-bit flags for the ith hit in list. All undefined bits are set
     * to zero. The following bits are defined:
     *
     * Bit Field  DAQ Meaning
     * 0-3        NCELL     (CMOS Cell address 0 to 15)
     * 4-5        FLAG1      Bit 4 = CGT ES16.
     *                       Bit 5 = CGT ES24.
     * 6-9        FLAG2      Bit 6 = Missed Count
     *                       Bit 7 = NC/CC Flag
     *                       Bit 8 = LGISELECT (1=long sample)
     *                       Bit 9 = CMOS ES16. */
    uint32_t pif;
    /* Uncalibrated times. */
    float pit;
    /* Uncalibrated high-gain, long int. charge. */
    float pihl;
    /* Uncalibrated high-gain, short int. charge. */
    float pihs;
    /* Uncalibrated log-gain, long int. charge. */
    float pilx;
    /* Uncalibrated unjittered time. */
    float pit0;
    /* PMT channel cell number. */
    uint32_t cell;
    /* CCC DAQ Circuit Number
     *     = 1024*card + 32*crate + channel.
     * Before uncalibration it holds the tube no. */
    uint32_t pin;
    /* TSLH of PMT. */
    uint32_t tslh;
    /* HCA information. */
    uint32_t hca;
    /* ECA validation status word. */
    uint32_t eca_val;
    /* PCA validation status word. */
    uint32_t pca_val;
    /* ANXX validation status word. */
    uint32_t anxx;
    /* ECA calibrated time (nsec). */
    uint32_t ept;
    /* ECA calibrated QHL (pedestal subtracted). */
    uint32_t ehl;
    /* ECA calibrated QHS (pedestal subtracted). */
    uint32_t ehs;
    /* ECA calibrated QLX (pedestal subtracted). */
    uint32_t elx;
    /* Non-walk corrected PMT time. */
    uint32_t pt1;
    /* Multiphoton PCA time. */
    uint32_t ptm;
    /* Multiphoton PCA PMT transit time RMS. */
    uint32_t ptms;
    /* Best charge (either QHS or QLX). */
    uint32_t qm;
    /* Best charge status word. */
    uint32_t qms;
    /* Charge correction for rate-dependent shifting. */
    uint32_t qrc;
} PMTBank;

void unpack_mcvx(uint32_t *data, MCVXBank *b);
void unpack_mctk(uint32_t *data, MCTKBank *b);
void unpack_ev(uint32_t *data, EVBank *b);
void unpack_pmt(uint32_t *data, PMTBank *b);

int isOrphan(aPmtEventRecord *pmtRecord);
void swap_int32(int32_t *val_pt, int count);
void swap_int16(int16_t *val_pt, int count);
int swap_PmtRecord(aPmtEventRecord *aPmtRecord, size_t size);
void swap_TrigRecord(struct TriggerInfo *aTrigRecord);
void swap_RunRecord(struct RunRecord *aRunRecord);

#endif