Calculations

CLast calculations are called to return an I3Particle as an early cascade fit, and to return I3CLastParams which details features of the Tensor of Inertia which are not encapsulated by the I3Particle.

The position of the I3Particle and the point about which the Tensor of Inertia is calculated is the Center of Gravity (CoG), called as I3Cuts::COG. The center of gravity is defined as

\(\vec{x} = \sum_i {a_i}^p \vec{r}_i\)

where \(\vec{r}_i\) is the position vector to a given DOM, and \(a_i\) is the amplitude of the DOM (the total charge) and \(p\) is the amplitude weight power which here is \(p = 1.0\).

The following calculations determine the direction via the Tensor of Inertia and also add how spherical the event is:

  • Tensor of Intertia (ToI) calculation:

    I3CLastCalculator::CalculateTOI takes the following variables: pulse series, IceCube geometry, and center of gravity (CoG). It then returns an I3Matrix containing the ToI.

    The formula for ToI is

    \(I^{ij} = \{\sum_k^N m_k (\delta^{ij}{R_k}^2 - {R_{k}}^{i} {R_{k}}^{j})\}/\sum_k^N m_k\)

    Where \(R_{ik}\) is the vector from the CoG to the kth mass, in our case the position of the kth DOM; and \({R_k}^2\) is the square of it’s length; \(m_k\) is the kth mass, in our case the total charge observed on the DOM. The calculation is for the \(I^{ij}\) component of the Tensor of Inertia which in our case is a 3x3 matrix.

    There are some variables which can modify the calculation: ampOpt_ and ampWeight_ default to 1.0, which results in predicable behaviour. If ampOpt_=0. then pulse charges between 0 and 2 are cast to 1.0. If ampWeight_ is set, each pulse is scaled to pow(charge,ampWeight_) (for example, ampWeight_=0. would cast all pulses regardless of charge to 1.0).

  • Diagonalising the ToI:

    I3CLastCalculator::diagonaliseToI is used as a means to access the eigenvalues of the tensor of intertia, and returns the eigenvector as a std::vector corresponding to the minimum eigenvalue. It uses the gsl library of matrix operations to perform this task.

    Note that the smallest eigenvalue will correspond to the long axis of the ‘object’ as this is the axis about which the ‘object’ could most freely rotate (i.e. has the least inertia).

    This calculation both dermines the eigenvector corresponding to the long axis, and also adds I3CLastParams to the frame which includes all 3 eigenvalues as well as the ratio evalratio \(e_1 / (e_1 + e_2 + e_3)\) which describes how spherical the event is.

  • Correct Direction

    I3CLastCalculator::CorrectDirection is used to remove the degeneracy (correct the calculation) of the two possible directions for the result eigenvector (forward or backward) by pointing the vector towards the direction where the OM hit time is latest from the CoG. To do this the covariance of the times projected along the vector \(l_i\) for DOM i is calculated weighted by the square of the charge. A negative covariance results in reversing the vector.

    The covariance calculation is done in two parts: a divisor \(C\) which is purely geometric:

    \(C = \sum_i {q_i}^2 \sum_i (l_i^2 {q_i}^2) - (\sum_i l_i {q_i}^2)^2\)

    and the result which is dependent on time \(t_i\):

    \(\sigma = [ \sum_i {q_i}^2 \sum_i (t_i l_i {q_i}^2) - \sum_i (l_i {q_i}^2) \sum_i(t_i {q_i}^2) ] / C\)

    If \(\sigma < 0\) then the eigenvector is reversed in order to point in the direction of later times.

In addition, to complete the seed requires a calculation of the time and energy of the event:

  • Time Calculation

    I3CLastCalculator::CalculateTimes is used to determine the time of the reconstruction particle where we currently have the CoG position and ToI eigenvector for direction.

    Basics of calculation: first we look at only DOMs within a radius rMax_, with a default of 300.0m get a “vertex time”. The vertex time is the time of an event at the CoG assuming the first hit observed on the given DOM was a direct photon from the CoG; i.e.

    \(t_{vertex} = t_{hit} - r/c_{ice}\)

    where \(r\) is the distance from the DOM to the CoG.

    Then, iterate over all other DOMs to check how many DOMs have a later vertex time than your current DOM (nDirect). To be a direct hit, the conditions are that the new hit trial vertex time \(t_{proj} = t_{hit2} - r_2 /c_{ice}\) must satisfy:

    \(t_{proj} > t_{vertex}\)

    and

    \(t_{proj} < t_{vertex} + t_{window}\).

    Here, \(t_{window}\) is the setting directHitWindow_ with a default of 200.0ns.

    To cover all bases, three times are stored: earliestThresholdTime for the earliest overall vertex time with \(n_{direct}\) times above the predefined threshold directHitThreshold_ with a default of 3 hits. (if satisfied, this will take precedence)

    timeOfMaxNDirect for the vertex time with the most \(n_{direct}\) times overall.

    earliestVertexTime as a fallback for the earliest vertex time regardless of threshold.

  • Energy Calculation

    I3CLastCalculator::CalculateEnergy_From_I3Hits and I3CLastCalculator::CalculateEnergy_From_AMHits are functions that covert the total charge observed to a corresponding energy given their own parameterisations.