.. _radcube_tutorial:

==============================
Tutorials on Using the Project
==============================

Below is a description on how to use this class to simulate an event. Note that there is an example event located in ``$I3_TESTDATA/radcube/100000``. This event was hand-made to some extent to be used for testing purposes. If you would like real simulated events, they are available on Cobalt (or you could make some yourself).

You can find example scripts to run in ``$I3_BUILD/radcube/resources/examples/`` which include calling typical module sequences as well as custom functions to augment the data for plotting. 

- The file ``ExampleScript_GetEFieldInVxB.py`` reads in a simulated event and converts the antenna locations to the :math:`\vec V \times \vec B` coordinate system, and makes a plot of the log of the peak amplitude of the Hilbert envelope.

- The file ``ExampleScript_SimulateRadioEvent.py`` is an example of how to perform a full simulation of an event and save it to a file. You will see that this is actually fairly similar thanks to the :ref:`RadioInjection` segment. 

There are two services which are important and unique to this project, ``I3AntennaResponse`` and ``I3ElectronicsResponse`` (see :ref:`services`). 

----------------------------
Use of Basic RadCube Modules
----------------------------

For a more detailed overview of what the various modules within :ref:`radcube`, see :ref:`modules`.

In general a simulation if you want to do it "by hand" would follow the code below

.. code-block:: python

   #First we import a library to handle python inputs smoothly
   import argparse
   parser = argparse.ArgumentParser()
   parser.add_argument('input', type=str, nargs='+', help='Input data files.')
   args = parser.parse_args()

   #import the standard libararies for icetray
   from icecube.icetray import I3Tray
   from icecube import icetray, dataio
   from icecube.icetray import I3Units

   #import radcube!
   from icecube import radcube

   tray = I3Tray()

   #We make one copy of all the services that we need
   #These will give you the "standard" setting for the respective services
   #The names of the services are returned so you can pass them to the respective module(s)
   antennaResponse = radcube.defaults.CreateDefaultAntennaResponse(tray)
   electronicsResponse = radcube.defaults.CreateDefaultElectronicsResponse(tray)
   randomService = radcube.defaults.CreateDefaultRandomService(tray)

   #Alternatively you could create them yourself with non-default parameters
   tray.AddService("I3ElectronicsResponseFactory", "MyNonDefaultER",
                    IncludeCables=True,
                    CableLength=100*I3Units.meter,
                    CableTemperature=274*I3Units.kelvin,
                    IncludeLNA=True,
                    IncludeRadioBoard=True,
                    IncludeTaxi=True,
                    AntennaType=antType,
                    IncludeElectronicNoise = False,
                    IncludeRadioBoardNoise = False,
                    InstallServiceAs="MyNonDefaultER"
                   )

   #Add the CoREAS reading module
   tray.Add("CoreasReader", "CoreasReader",
            DirectoryList=args.input,      #List of CoREAS output directories (one per shower))
            MakeGCDFrames=True,            #Make a GCD frames from simulations
            MakeDAQFrames=True,            #Make the Q frames from simulations
           )

   #Inject the electric fields into the antennas
   tray.AddModule("ChannelInjector", "_ChannelInjector",
                  InputName="CoREASEFieldMap",  #This is the default name that the CoreasReader puts in the frame
                  OutputName="RawVoltageMap",
                  AntennaResponseName=antennaResponse
                 )

   #Convolve with the electronics response
   #Here we will use the default response we created above, but you can also use
   #the non-default one as well. Plug and play!
   tray.AddModule("ElectronicResponseAdder", "AddElectronics",
                  InputName="RawVoltageMap",
                  OutputName="FoldedVoltageMap",
                  ElectronicsResponse=electronicsResponse
                 )

   #Make yourself a P Frame (CoreasReader only makes GCD and Qs)
   tray.AddModule("I3NullSplitter","splitter",
                  SubEventStreamName="RadioEvent"
                 )

   #Filter the voltages to only look at a narrower band in frequency space
   #Here we will filter between 110 and 190 MHz using the Butterworth algorithm
   tray.AddModule("BandpassFilter", "BandpassFilterButter",
                  InputName="RadioEvent",
                  OutputName="Filtered",
                  FilterType=radcube.eButterworth,
                  FilterLimits=[110*I3Units.megahertz,190*I3Units.megahertz],
                  ButterworthOrder=radcube.constants.butterworthOrder
                 )

   #Add any other modules you want here for: analyzing, plotting, saving, etc.

   tray.Execute()

Likewise one would run this code as

.. code-block:: console

    $ python the-script-defined-above.py $I3_TESTDATA/radcube/100000

----------------------------------------
Full Radio Reconstruction Using Segments
----------------------------------------

In general, there are two main segments for :ref:`radcube`, the :ref:`RadioInjection` and :ref:`RadioReconPrep` (defined in ``./radcube/python/segments/``). The first reads in coreas simulations, creates GCD frames, as needed, and makes one Q frame per passed in event. The second takes simulated files (expected to be in a format that comes from the Pole) and prepares them for reconstruction. This includes removing various antennas from the event, estimating the shower properties, etc.

These two segments will be kept up to date such that an "official" (or at least recommended) simulation would instead be the more simple:

.. code-block:: python

    import argparse
    import os

    parser = argparse.ArgumentParser()
    parser.add_argument('--output', dest='outfile', type=str, default=str(os.getenv('I3_BUILD'))+'/simulated-radio-event.i3.gz', help='Select to add radio to input gcd file')
    parser.add_argument('input', type=str, nargs='+', help='Input data files.')
    args = parser.parse_args()

    from icecube.icetray import I3Tray
    from icecube import icetray, radcube, dataio
    from icecube.icetray import I3Units
    from icecube.icetray.i3logging import *
    from icecube.radcube.defaults import *

    tray = I3Tray()

    #These lines create basic response services. See radcube.defaults for more
    #The functions place them in the frame and give you back the name
    antennaResponse = CreateDefaultAntennaResponse(tray)
    electronicsResponse = CreateDefaultElectronicsResponse(tray)
    randomService = CreateDefaultRandomService(tray)

    #This segment reads in the CoREAS event, adds noise, convolves it with the
    #antenna/electronics response, "digitizes" it, and puts it in the frame.
    tray.AddSegment(radcube.segments.RadioInjection, 'RadioInjection',
                     InputFiles = args.input,
                     AntennaServiceName = antennaResponse,
                     ElectronicServiceName = electronicsResponse,
                     RandomServiceName = randomService,
                     OutputName = "InjectedVoltageMap"
                    )

    #Make yourself a P Frame (CoreasReader only makes GCD and Qs)
    tray.AddModule("I3NullSplitter","splitter",
                  SubEventStreamName="RadioEvent"
                 )

    #This segment takes the simulated data, converts from ADC bins to a 
    #voltage, removes electronics effects, and estimates shower parameters
    tray.AddSegment(radcube.segments.RadioReconPrep, "ReconPrep",
                    ElectronicServiceName = "electronicResponseName",
                    InputName = "InjectedVoltageMap"
                   )

    tray.AddModule("I3Writer","i3writer",
                   filename=args.outfile,
                   # streams=[icetray.I3Frame.DAQ,icetray.I3Frame.Physics]
                   # streams=[icetray.I3Frame.Geometry,icetray.I3Frame.Calibration, icetray.I3Frame.DetectorStatus]
                   )

    tray.Execute()

Again, this sctipt can be run by passing in an event

.. code-block:: console

   $ python the-script-defined-above.py $I3_TESTDATA/radcube/100000

You can also pass in an output location to put the file

.. code-block:: console

   $ python the-script-defined-above.py --output ./myfile.i3.gz $I3_TESTDATA/radcube/100000