Usage

IceTray Service

The intended usage of rpdf is to provide the likelihood service for a Gulliver reconstruction. As shown in the following example:

tray.AddService("I3RecoLLHFactory", "spe_likelihood",
   InputReadout="InIcePulses",
   JitterTime=15*I3Units.ns,
   NoiseProbability=10.*I3Units.hertz
)

tray.Add("I3SimpleFitter",
  OutputName="SPESingleFit",
  SeedService="seed",
  Parametrization="simpletrack",
  LogLikelihood="spe_likelhood",
  Minimizer="minuit",
)

In the example an instance of I3RecoLLH with the name "spe_likelihood" is placed in the tray’s context by I3RecoLLHFactory. This example reads an I3RecoPulsesSeriesMap from the frame and uses a jitter time of 15 ns and a noise rate of 10 Hz. A Gulliver reconstruction I3SimpleFitter is then added to the tray and passed the name of the rpdf service to the LogLikelihood parameter.

Python

The low level functions are also available via a python interface. An example calculation is below:

from icecube import rpdf

icemodel = rpdf.H2
peprob = rpdf.FastConvolutedPandel()
jitter = 15*I3Units.ns
noise =  10*I3Units.hz
llhtot = 0

for omkey, pulseseries in self.pulses:

  #find the position of the OM
  ompos = geometry[omkey].position

  #calculate the geometrical parameters tgeo and deff
  geo_params = rdf.muon_geometry(ompos,part,icemodel)

  #residual time is the pulse time minus tgeo
  t_res= pulseseries[0].time-geo_params.first

  #calculate the convoluted Pandel function of t_res
  llh = peprob.pdf(t_res,geo_params.second,jitter,icemodel)

  #add noise and sum the log of the likelihoods
  llhtot += math.log(self.noise + llh)

print llhtot

This is equivalent to calculating the SPE1st likelihood in I3RecoLLH.

A easier method if you want to calculate the likelihood of a single event hypothesis is to instantiate a I3RecoLLH object in python land:

recollh = rpdf.I3RecoLLH("SRTHVInIcePulses","SPE1st","GaussConvoluted",
                         15*I3Units.ns,10.*I3Units.hertz,rpdf.H2)

particle = phys["SPEFitSingle_TWHV"]
hypothesis = gulliver.I3EventHypothesis(particle)
recollh.set_geometry(i3geometry)
recollh.set_event(physics_frame)