aboutsummaryrefslogtreecommitdiff
path: root/src
AgeCommit message (Collapse)Author
2018-10-12skip PMTs which weren't hit for the fast likelihood calculationtlatorre
2018-10-06pre { line-height: 125%; } td.linenos .normal { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; } span.linenos { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; } td.linenos .special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; } span.linenos.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; } .highlight .hll { background-color: #ffffcc } .highlight .c { color: #888888 } /* Comment */ .highlight .err { color: #a61717; background-color: #e3d2d2 } /* Error */ .highlight .k { color: #008800; font-weight: bold } /* Keyword */ .highlight .ch { color: #888888 } /* Comment.Hashbang */ .highlight .cm { color: #888888 } /* Comment.Multiline */ .highlight .cp { color: #cc0000; font-weight: bold } /* Comment.Preproc */ .highlight .cpf { color: #888888 } /* Comment.PreprocFile */ .highlight .c1 { color: #888888 } /* Comment.Single */ .highlight .cs { color: #cc0000; font-weight: bold; background-color: #fff0f0 } /* Comment.Special */ .highlight .gd { color: #000000; background-color: #ffdddd } /* Generic.Deleted */ .highlight .ge { font-style: italic } /* Generic.Emph */ .highlight .ges { 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 */
#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);
}
at only works when the x array is uniform. Since the wavelengths are not spaced uniformly, we have to use the GSL interpolation routines. 2018-09-26speed up fast likelihood calculationtlatorre This commit updates the fast likelihood calculation to use the identity sin(a-b) = sin(a)*cos(b) - cos(a)*sin(b) to speed up the fast likelihood calculation. 2018-09-26speed up fast likelihood calculationtlatorre This commit speeds up the fast likelihood calculation by avoiding calls to trigonometric functions where possible. Specifically we calculate sin(a) = sqrt(1-pow(cos(a),2)); instead of sin(a) = sin(acos(cos(a))); 2018-09-26increase number of points in cos(theta) interpolation to 1000tlatorre 2018-09-25update indirect scattering PDF start timetlatorre Currently the PDF for scattered light is modelled as a flat distribution starting at some time t. Previously I was using the mean hit time for all PMTs, however this should really be a flat distribution in the time *residual* after the main peak. Therefore, the PDF now starts at the estimated time for direct photons. 2018-09-25update likelihood calculation to use PMT_TTS macrotlatorre I accidentally hardcoded the single PE TTS to 1.5 ns in the likelihood calculation. 2018-09-25update integration bounds in likelihood calculationtlatorre This commit updates the bounds of the track integration in the likelihood function to integrate up to 1 meter around the point at which the PMT is at the Cerenkov angle from the track. This fixes an issue I was seeing where a *very* small change in the fit paramters would cause the likelihood to jump by a large amount. I eventually tracked it down to the same issue I was seeing before which I solved by splitting up the integration into two intervals. However that fix did not seem to completely fix the issue. Based on initial tests with 500 MeV muons, this fix seems to do a much better job. 2018-09-25increase maxeval to 20 for the "quick" minimization phasetlatorre When testing out the fitter on 500 MeV muons, there was at least one event which started out at a position very far from the true position. This event had a secondary electron like ring which is what I suspect caused the fit to start out in a position far from the true position. This fix correctly starts the minimization close to the true position. In the future I should look at updating get_direction() so that it finds the largest ring direction instead of just doing a weighted average of all the vectors from the position to the PMTs. 2018-09-21update likelihood function to include the probability of the path coefficientstlatorre 2018-09-21split up the track integral into two intervalstlatorre This commit updates the likelihood calculation to split up the track integral into two intervals in some cases. I noticed when fitting some events that the likelihood value would change drastically for a very small change in the fit parameters. I eventually tracked it down to the fact that the track integral was occasionally returning a very small charge for a PMT which should have a very high charge. This was happening because the region of the track which was hitting the PMT was very small and the cquad integration routine was completely skipping it. The solution to this problem is a bit of a hack, but it seems to work. I first calculate where along the track (for a straight track) the PMT would be at the Cerenkov angle from the track. If this point is somewhere along the track then we split up the integral into two intervals: one going from the start of the track to this point and the other from the point to the end of the track. Since the cquad routine always samples points near the end of the intervals this should prevent it from completely skipping over the point in the track where the integrand is non-zero. 2018-09-21increase MAX_PE to 10000tlatorre 2018-09-20delete unused variable distancetlatorre 2018-09-20don't include the OWL PMTs in the likelihood calculationtlatorre For some reason the OWL tubes have 9999.00 for the x, y, and z coordinates of the normal vector in the PMT file. For now, I'm just going to remove them from the likelihood calculation. 2018-09-20make sure direction vector is normalized in path_eval()tlatorre 2018-09-20add particle id code to output filetlatorre 2018-09-20add git SHA1 hash to output filetlatorre 2018-09-20add time elapsed to the output filetlatorre 2018-09-20add a command line option to only fit the first eventtlatorre This commit adds the ability to run the fit program with the --skip-second-event flag to only fit the first event after a MAST bank. This way we avoid fitting secondaries like Michel electrons when fitting MC events. 2018-09-20add absorption lengths for D2O and H2O weighted by the Cerenkov spectrum and ↵tlatorre the quantum efficiency 2018-09-19change output file format to YAMLtlatorre 2018-09-18stop fitting when the likelihood difference is less than 1e-2tlatorre 2018-09-18speed likelihood calculation up a bittlatorre 2018-09-18free muon_energy struct in test_muon_get_energy()tlatorre 2018-09-18fix heap overflow in interp1d()tlatorre This commit fixes a potential heap overflow in interp1d() which could occur if x was very close to the last value in the xp array. The bounds check is now performed on the index rather than the x values. 2018-09-18fix typotlatorre 2018-09-18update PMT TTS and dark rate from pmt_response.dattlatorre 2018-09-18fix typotlatorre 2018-09-18fix PMTR bank offsetstlatorre 2018-09-18fix offsettlatorre 2018-09-18add functions to get the PMT reflectivity from the PMTR banktlatorre 2018-09-18add free_charge() to free memory used to interpolate the charge distributionstlatorre