Photon Simulation Tutorial » History » Version 3

« Previous - Version 3/46 (diff) - Next » - Current version
Jonathan Insler, 01/05/2016 10:14 AM
corrected 'dunecode' to 'dunetpc'

Photon Simulation Tutorials

This tutorial will explain how to use the basic tools for photon simulation. It assumes a basic working knowledge of setting up and using LArSoft (e.g. have used the AnalysisExample). You will learn:

  • How to simulate and examine the photon signals on photon detectors in a physics events.
  • (How to examine an existing photon library.)
  • (How to generate a new photon library.)

Tutorial 0: Setup your environment

Create a test release with the latest version of LArSoft. You will need to check out (mrb g) the following packages from develop:

  • dunetpc

Set up your environment, build and install the packages. (Not sure how? Look here or here)

Tutorial 1: Generate photon signals

Go to this directory in the dunetpc package which contains the tutorial files:


Take a look a the file dune35t_optical_tutorial_sim.fcl. It is a basic prodsingle example that will produce photon signals from a muon and put the information into root trees we can examine. There are several key pieces required to make the photon simulation work. The first is to include the optical detector fhicl files:

#include "photpropservices_dune.fcl" 
#include "opticaldetectormodules_dune.fcl" 

The next is to include pmtresponse in the analyzers and add it to the analyzeIt end path:

   pmtresponse: @local::dune35t_simphotoncounter
   analyzeIt:  [ AnalysisExample, pmtresponse]
 end_paths:     [analyzeIt, stream1]  

Then the photon library must be set up using the default 35 ton configuration:

# DUNE 35ton w/o wires and mesh
services.user.PhotonVisibilityService: @local::dune35t_photonvisibilityservice

Finally, the optical physics processes need to be turned on in Geant4:

# set quantum efficiency supressed scint yield to 0.03 * 24000
services.user.LArProperties.ScintYield: 24000

# enable optical physics in LArG4
services.user.LArG4Parameters.EnabledPhysics: [ "Em",
                                                "NeutronTrackingCut" ]

# enable this custom physics list
services.user.LArG4Parameters.UseCustomPhysics: true

# disable cerenkov light
services.user.LArProperties.EnableCerenkovLight: false

All these pieces are set up ahead of time in the sim fhicl file. To begin, just try simulating an event:

lar -c dune35t_optical_tutorial_sim.fcl -n 1

This will produce two files:


The first contains the ART event records and the second contains simple root trees we can examine. Open up the second file and cd into the pmtresponse directory which contains the trees generated by the pmtresponse modules.

Attaching file dune35t_optical_tutorial_sim_hist.root as _file0...
root [1] .ls
TFile**         dune35t_optical_tutorial_sim_hist.root
 TFile*         dune35t_optical_tutorial_sim_hist.root
  KEY: TDirectoryFile   largeant;1      largeant (LArG4) folder
  KEY: TDirectoryFile   AnalysisExample;1       AnalysisExample (AnalysisExample) folder
  KEY: TDirectoryFile   pmtresponse;1   pmtresponse (SimPhotonCounter) folder
root [2] pmtresponse->cd()
root [3] .ls
TDirectoryFile*         pmtresponse     pmtresponse (SimPhotonCounter) folder
 KEY: TTree     AllPhotons;1    AllPhotons
 KEY: TTree     DetectedPhotons;1       DetectedPhotons
 KEY: TTree     OpDets;1        OpDets
 KEY: TTree     OpDetEvents;1   OpDetEvents

These trees contain the true photon information from the muon we just simulated. AllPhotons contains one entry per simulated photon per event which reaches a PMT, storing the event in which the photon was produced, the optical channel, and the time of the photon. In our simulation the response of each photon detector is assumed to be perfect, so the DetectedPhotons tree contains the same information. OpDets has one entry per photon detector per event containing the optical detector channel and the number of photons that reached it.

To start, make some simple distributions from these trees. For example:
  • What is the photon time distribution?
  • Which optical channel had the largest photon signal?

Tutorial 2: Compare photon signals

In this tutorial, we will generate muons at range of different positions in the detector and use the changing signals in each of the photon detectors to determine the location of each.

Here we will use textfilegen instead of singlesgen as the generator. This allows a repeatable set of input particles to be used. The text files which will serve as inputs can be found in the vectors subdirectory and they are called by the three scan fhicl files:


Each fhicl file/vector combination will produce 20 muons at a range of positions in x, y, and z, respectively. Let's start by scanning in the y direction.

lar -c dune35t_optical_tutorial_yscan.fcl

That will run for about 4 minutes -- get a cup of coffee or check your e-mail. When it's done, you'll again have 2 root files:


These have the same elements as in Tutorial 1. Now, however, there are more events and so more entries in all the trees. We want to try to extract the information from those trees in order to make sense of the simulated data. Included in the tutorial directory is a ROOT macro that shows you how to plot the signals from a single photon detector vs. the true starting position of the muons. The photon detector signal comes from the pmtresponse/OpDets tree while the true muon starting position comes from the AnalysisExample/AnalysisExampleSimulation tree.

Run it:

root -l Tutorial2.C

You should get a plot that looks like this:

From this plot, you can probably guess where (top, middle or bottom) this photon detector is placed.

You can save the plot you just made to a file by typing the following at the ROOT command prompt:

Now, try the following:
  • Take a look at the other Photon Detectors. See if you can determine the location of each of them.
  • Once you've located each of the PDs in y, try doing the scans in x and z. Again, try to locate the PDs in each dimension.

Tutorial 3: How to examine a photon library directly

The photon library store the pre-recorded sensitivity of each photon detector to each voxel (3-dimensional pixel) of the detector. While the library is slow to produce, it then makes running photon simulation fast. The contents of the library can be examined directly, which we demonstrated in this tutorial.

We will be analyzing the default 35ton photon library already in the repository, the same one you used in the first two tutorials. The library is located in the larsoft_data package. You likely don't have it checked out so to look at it directly you'll need to go to the base product (If this path doesn't work, look at your $FW_SEARCH_PATH variable and look for larsoft_data.):


The parent directory there has a README listing the properties of all the committed libraries. If you open up that file in root, you will see a single tree called PhotonLibraryData. Let's take a look at the contents of that tree, just printing out the values in the first entry:

Attaching file /grid/fermiapp/products/larsoft/larsoft_data/v0_02_00/PhotonPropagation/LibraryData/lib11225235.root as _file0...
root [1] .ls
TFile**         /grid/fermiapp/products/larsoft/larsoft_data/v0_02_00/PhotonPropagation/LibraryData/lib11225235.root
 TFile*         /grid/fermiapp/products/larsoft/larsoft_data/v0_02_00/PhotonPropagation/LibraryData/lib11225235.root
  KEY: TTree    PhotonLibraryData;1     PhotonLibraryData
root [2] PhotonLibraryData->Show(0)
======> EVENT:0
 Voxel           = 0
 OpChannel       = 7
 Visibility      = 0.000233333

(This is actually from an older tree so your numbers will vary.)

We see that this entry has a Voxel, and optical channel, and a Visibility, so it tells us the sensitivity of this channel to this voxel. You can poke around a few more channels, but this tree has a lot of entries (> 5 million) so this is a pretty inefficient way to learn anything. Instead let's take advantage of an existing LArSoft module just for this purpose: PhotonLibraryAnalyzer.

The relevant fhicl file here is dune35t_optical_tutorial_libanalysis.fcl. Run it and let's look at the output.

lar -c dune35t_optical_tutorial_libanalysis.fcl

This produces a file called PhotonLibraryFile_dune35texample_hists.root. It contains projections of the photon library in x, y, and z (summed across the axis not plotted) as well as numerous slices through the unplotted axis. Start with the summed projection in the x direction. This should give you a clear picture of the layout of the photon detectors. Now, step through a few of the slices to see how the sensitivity in y and z changes as a function of x. Sample root commands below:

Attaching file PhotonLibraryFile_dune35texample_hists.root as _file0...
root [1] .ls
TFile**        PhotonLibraryFile_dune35texample_hists.root
 TFile*        PhotonLibraryFile_dune35texample_hists.root
  KEY: TDirectoryFile    libana;1    libana (PhotonLibraryAnalyzer) folder
root [2] libana->cd()
root [3] .ls
root [4] XProjection->Draw("colz")
root [5] projX5->Draw("colz")
root [6] projX10->Draw("colz")
root [7] projX20->Draw("colz")
root [8] projX40->Draw("colz")

Once you've done that, explore the other projections a bit to get a fuller sense of the 35 ton geometry as it is seen by the photon detectors.

Tutorial 4: How to generate digitized photon detector signals

In this tutorial, we will take simulated photons, like those produced in the first and second tutorial, and pass them through a simulation of the SiPM and front-end electronics which "digitize" the output to produce signals which look more like output a real detector. The fhicl file which simulates muons and then digitizes the photon signals is called dune35t_optical_tutorial_digi.fcl and can be found in dunetpc/dune/PhotonPropagation/Tutorial/.

The digitization is performed by the OpMCDigi module, which is called by this line in the producers section of the fhicl file:

opdigi:    @local::dunefd_opdigi

We can examine the output of the digitizer with a pre-built analyzer module, called from the analyzers section:

opdigiana:   @local::dune35t_opdigiana

Now, run the fhicl file and let's look at the output. This fhicl file uses the output of the first tutorial as an input, which you specify with the -s <name.root> command. If you renamed the output file from the first tutorial, use that name instead of the name listed below.

lar -c dune35t_optical_tutorial_digi.fcl -s dune35t_optical_tutorial_sim_gen.root

As usual, it will produce two files, dune35t_optical_tutorial_digi_hist.root, and the same with gen instead of hist. Open up the hist file in root, and let's look at the contents of the opdigiana directory:

root [0]
Attaching file dune35t_optical_tutorial_digi_hist.root as _file0...
root [1] .ls
TFile**         dune35t_optical_tutorial_digi_hist.root
 TFile*         dune35t_optical_tutorial_digi_hist.root
  KEY: TDirectoryFile   daq;1   daq (SimWireDUNE35t) folder
  KEY: TDirectoryFile   largeant;1      largeant (LArG4) folder
  KEY: TDirectoryFile   AnalysisExample;1       AnalysisExample (AnalysisExample) folder
  KEY: TDirectoryFile   pmtresponse;1   pmtresponse (SimPhotonCounter) folder
  KEY: TDirectoryFile   opdigiana;1     opdigiana (OpDigiAna) folder
root [2] opdigiana->cd()
root [3] .ls
TDirectoryFile*         opdigiana       opdigiana (OpDigiAna) folder
 KEY: TH1D      event_1_opdet_0;1
 KEY: TH1D      event_1_opdet_1;1
 KEY: TH1D      event_1_opdet_2;1
 KEY: TH1D      event_1_opdet_12;1
 KEY: TH1D      event_1_opdet_13;1
 KEY: TH1D      event_1_opdet_14;1
 KEY: TH1D      event_1_opdet_15;1
 KEY: TH1D    event_1_opdet_95;1

You'll see one histogram per OpDet channel per event (channels with no data are skipped). Let's take a look at one of the channels.

root [4] event_1_opdet_12->Draw("hist")

You should see a charge vs. time distribution looking something like this, though you probably will need to zoom in on the x-axis in order to see the photon signal which will be in the first few us close to 0.

Tutorial 5: How to run photon detector reconstruction and examine the output

In this tutorial, we will take the digitized photon signals produced in Tutorial 4 and run the photon detector reconstruction algorithms on them. The goal of the reconstruction is effectively to work backwards to approximate the truth based on the observed signals. The first step goes from the digitized waveforms to the individual or groups of photons arriving at the photon detectors at a particular time, called OpHits. The second step groups together simultaneous OpHits across multiple detectors into OpFlashes, which in principle correspond to a source of scintillation light at a particular place and time in the detector, such as an electron.

The fhicl file which simulates muons and then digitizes the photon signals is called dune35t_optical_tutorial_reco.fcl and can be found in dunetpc/dune/PhotonPropagation/Tutorial/. This tutorial runs on the output of Tutorial 4, so please do that tutorial first.

The reconstruction is performed by the OpFlash module which builds both the OpHits and OpFlashes:

# photon detector reconstruction
  opflash:            @local::dune35t_opflash

We can examine the reconstructed objects with a pre-built analyzer module, called from the analyzers section:

opflashana:   @local::dune35t_opflashana

Now, run the fhicl file and let's look at the output. Note that this fhicl file runs over the output of the previous tutorial, so if you renamed the output file from that tutorial specify the new source name with the -s <new name.root> command.

lar -c dune35t_optical_tutorial_reco.fcl -s dune35t_optical_tutorial_digi_gen.root

As usual, it will produce two files, dune35t_optical_tutorial_reco_hist.root, and the same with gen instead of hist. Open up the hist file in root, and let's look at the contents of the opflashana directory:

root [0]
Attaching file dune35t_optical_tutorial_reco_hist.root as _file0...
root [1] .ls
TFile**        dune35t_optical_tutorial_reco_hist.root
 TFile*        dune35t_optical_tutorial_reco_hist.root
  KEY: TDirectoryFile    opflashana;1    opflashana (OpFlashAna) folder
  KEY: TDirectoryFile    AnalysisExample;1    AnalysisExample (AnalysisExample) folder
  KEY: TDirectoryFile    pmtresponse;1    pmtresponse (SimPhotonCounter) folder
  KEY: TDirectoryFile    opdigiana;1    opdigiana (OpDetDigiAnaDUNE) folder
root [2] opflashana->cd()
root [3] .ls
TDirectoryFile*        opflashana    opflashana (OpFlashAna) folder
 KEY: TTree    FlashBreakdownTree;1    FlashBreakdownTree
 KEY: TTree    PerOpHitTree;1    PerOpHitTree
 KEY: TTree    PerFlashTree;1    PerFlashTree
 KEY: TTree    FlashHitMatchTree;1    FlashHitMatchTree

Four root trees are made which give you different ways to examine the output of reconstruction.


PerOpHitTree contains one entry per OpHit across all optical detectors and all events. Take a look at the information available for each OpHit using the Show method:

root [4] PerOpHitTree->Show(0)
======> EVENT:0
 EventID         = 1
 HitID           = 0
 OpChannel       = 0
 PeakTimeAbs     = 0.226562
 PeakTime        = 0.226562
 Frame           = 0
 Width           = 3.48438
 Area            = 31376
 Amplitude       = 8
 PE              = 23.591
 FastToTotal     = 0
The most important information for each hit is contained in these four variables:
  • EventID: The event this hit occurred in. In this tutorial we have only produced one event, so all entries will have EventID 1.
  • OpChannel: The optical channel this hit was detected on. While there are only 8 optical detectors in the 35ton prototype we have simulated, each optical detector is read out by multiple optical channels so this will range from 0 to 96.
  • PeakTime: The time the hit occurred.
  • PE: The amount of light, in units of photo electrons, observed in this hit.

Try drawing the PE and PeakTime distributions. Examples are included below, though yours may look a little different due to statistical fluctuations.

The PE distribution is peaked towards zero since smaller hits are more common than larger hits. It is spiky both because the statistics are low and because the SiPM detectors are sensitive to single photons, meaning one, two, or three photon hits can be distinguished from one another, which will appear as separate peaks on the PE spectrum. You can try simulating more events and drawing this distribution with more bins to see this structure.

The time distribution will have two components. One set of hits that appear promptly near 0 (because this is simulation) and a scatter of hits at later times. This is a property of scintillation light in argon, which has these two (prompt and late) components.


PerFlashTree contains information for each reconstructed flash. Since we produced only one event with only one particle, this tree will likely only have one or two entries in it. Let's take a look at one of those entries:

======> EVENT:0
 EventID         = 1
 FlashID         = 0
 YCenter         = 21.5174
 ZCenter         = 90.4797
 YWidth          = 61.4153
 ZWidth          = 135.424
 FlashTime       = 0.266388
 AbsTime         = 0.266388
 FlashFrame      = 1
 InBeamFrame     = 0
 OnBeamTime      = 1
 TotalPE         = 321.573

Like the OpHit objects, the OpFlash has a time, FlashTime, and a total amount of light in TotalPE. Since it is built up from multiple detectors, it doesn't have a single channel but instead has a reconstructed position (and width) in Y and Z made up from the weighted averages of the positions of the photon detectors.

Other Trees

FlashHitMatchTree contains one entry for each hit which is a part of a flash is used to dig into the hits which make up a particular flash. It includes some information from the PerOpHit and PerFlash trees, including the FlashID and HitID which allow the relevant entries in the other trees to be identified.

FlashBreakdownTree contains one entry for each optical detector in each flash. It contains all the information in the PerFlash tree, plus an OpChannel and an NPe variable which is amount of light seen on that particular optical channel in that flash.

Tutorial 6: How to examine and reformat raw photon detector data

First, you need to obtain some raw data to work with. You can get a raw data file from vertical slice (warm), reading out LED flashes from the SAM database (need to be on the Fermilab computers):

ifdh_fetch -e dune v35t_r0000101_s0001_test_raw.root

Then, the first way to look at the data is to dump information out to the screen. The fhicl file to do this can be found in dunetpc/dune/daqinput35t:

lar -c SSPEventDump.fcl -s v35t_r0000101_s0001_test_raw.root

Noting we use the "-s" command to specify the source root file. This will dump information about each bit of entry of SSP data to the screen, like below. This is from one of the entries that has some data in it, so you can see some example ADC values picked out of the recorded waveform on the Data Values line.

Begin processing the 18th record. run: 101 subRun: 1 event: 18 at 15-Jan-2015 14:34:13 CST

Run 101, subrun 1, event 18 has 1 fragment(s) of type PHOTON

SSP fragment 0 has total size: 49682 and run number: 999 with 49680 total ADC values

Header:                             2655784300
Length:                             640
Trigger type:                       0
Status flags:                       0
Header type:                        0
Trigger ID:                         58924
Module ID:                          2543
Channel ID:                         3
Sync delay:                         640
Sync count:                         12426
Peak sum:                           12
Peak time:                          0
Prerise:                            49682
Integrated sum:                     255744
Baseline:                           0
CFD Timestamp interpolation points: 43690 43690 1035 0
Internal interpolation point:       0
Internal timestamp:                 70368705118224

Data values: 43690 43690 1035 0 0 16 64940 16383 0 0 12904 61184 57573 44553 2745 5175 65517 65534 12 33 ...

Event ADC average is (from counter):   4.68166e+07
Event ADC average is (from histogram): 681.228

DAQ......... | daq......... | TOY1................. | std::vector<artdaq::Fragment> | artdaq::Fragments.... | ...?
DAQ......... | daq......... | PHOTON............... | std::vector<artdaq::Fragment> | artdaq::Fragments.... | ...1
DAQ......... | daq......... | TPC.................. | std::vector<artdaq::Fragment> | artdaq::Fragments.... | ...?
DAQ......... | daq......... | missed............... | std::vector<artdaq::Fragment> | artdaq::Fragments.... | ...?
DAQ......... | daq......... | unidentified......... | std::vector<artdaq::Fragment> | artdaq::Fragments.... | ...?
DAQ......... | daq......... | TOY2................. | std::vector<artdaq::Fragment> | artdaq::Fragments.... | ...?
DAQ......... | daq......... | TRIGGER.............. | std::vector<artdaq::Fragment> | artdaq::Fragments.... | ...?

Total products (present, not present): 7 (1, 6).

Next, we will reformat the data into an ART data product, OpDetPulse, that we can use with existing larsoft tools like OpDigiAna. The relevant fhicl files are in the same directory as above, dunetpc/dune/daqinput35t.

lar -c RunSSPToOffline.fcl -s v35t_r0000101_s0001_test_raw.root -o SSPEventToOffline_v35t_r0000101_s0001_test_raw.root

This will produce a root file whose filename you defined with the -o argument. This file will contain an art data product with the waveforms extracted from the raw SSP data. Open it up in a simple root session, and create a TBrowser:

root[1] new TBroswer()

We cannot access all of the information, but we can at least see that the data product is there. Open up the root file, and then click down into the Events tree. Inside, you will see a branch called:


This is the data product containing our photon detector data. To take a look at the waveforms we can use the same OpDigiAna module we used to look at the output of OpMCDigi since they both produce the same data product. However, to read in the raw data, we need to tell OpDigiAna the correct name of the data product to read in. This is done by changing two settings in the analysis module:

physics.analyzers.opdigiana.InputModule: ssptooffline
physics.analyzers.opdigiana.InstanceName: offlinePhoton

A fhicl file with these settings already in place can be found in the tutorial directory. Run it, setting the input directory and the TFileService (histogram) output file.

lar -c dune35t_optical_tutorial_rawdata_analysis.fcl -s  SSPEventToOffline_v35t_r0000101_s0001_test_raw.root -T opdigiana_hist.root

You can then look at the output in opdigiana_hist.root and see the waveforms from the vertical slice, which should look something like this:

Bonus Tutorial: How to generate your own photon library

This tutorial will not be necessary for everyone, so you probably don't need to go through it until you need it. The hope is that libraries will only need to be remade occasionally.