Project

General

Profile

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 working ara with the latest version of LArSoft and clone dunetpc. Set up your environment, build and install the packages. (This procedure not look familiar? Start here or here). Note that we use a specific version of larsoft known to work with these tutorials.

ups list -aK+ larsoft
setup larsoft v<latest version> -qdebug:e17
mkdir testRel_tutorial
cd testRel_tutorial
mrb newDev
source localProducts*/setup
cd srcs
mrb g dunetpc
mrbsetenv
mrb i

Tutorial 1: Generate photon signals

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

dunetpc/dune/PhotonPropagation/Tutorial

Explanation of the configuration

Take a look a the file dune1x2x6_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" 

Then we choose to use the 1x2x6 workspace geometry, and set up the corresponding "PhotonVisibilityService" which handles the photon library.

services.Geometry:                @local::dune10kt_1x2x6_geo
services.PhotonVisibilityService: @local::dune10kt_1x2x6_photonvisibilityservice

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

 analyzers:
 {
   pmtresponse: @local::dunefd_simphotoncounter
...
   analyzeIt:  [ AnalysisExample, pmtresponse]
...
 end_paths:     [analyzeIt, stream1]  
}

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

# Custom Physics List
services.LArG4Parameters.UseCustomPhysics: true
services.LArG4Parameters.EnabledPhysics: [ "Em",
                                           "FastOptical",
                                           "SynchrotronAndGN",
                                           "Ion",
                                           "Hadron",
                                           "Decay",
                                           "HadronElastic",
                                           "Stopping",
                                           "NeutronTrackingCut" ]

All these pieces are set up ahead of time in the sim fhicl file.

Running the tutorial

To begin, just try simulating an event:

lar -c dune1x2x6_optical_tutorial_sim.fcl -n 1

Some things you might notice in the output:

%MSG-i GeometryCore:  Early 07-Sep-2017 10:24:06 CDT JobSetup
New detector geometry loaded from
        /dune/app/users/ahimmel/testrels/testRel_tutorial/localProducts_larsoft_v06_48_00_e14_prof/dunetpc/v06_48_00/gdml/dune10kt_v1_1x2x6.gdml
        /dune/app/users/ahimmel/testrels/testRel_tutorial/localProducts_larsoft_v06_48_00_e14_prof/dunetpc/v06_48_00/gdml/dune10kt_v1_1x2x6_nowires.gdml

This is normal. The photon simulation occurs in a "parallel world" where the TPC wires have been removed to speed up the simulation.

Geant4 simulated 8 MC particles, we keep 8 .
%MSG
Begin processing the 5th record. run: 1 subRun: 0 event: 5 at 12-Sep-2017 09:25:26 CDT

This warning is harmless. It simply means there are areas where particles can pass where the electrons are not recorded, such as outside the detector or in the gaps between the drift regions.

Examining the output

This will produce two files:

dune1x2x6_optical_tutorial_sim_gen.root
dune1x2x6_optical_tutorial_sim_hist.root

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.

$ root dune1x2x6_optical_tutorial_sim_hist.root
   ------------------------------------------------------------
  | Welcome to ROOT 6.08/06                http://root.cern.ch |
  |                               (c) 1995-2017, The ROOT Team |
  | Built for linuxx8664gcc                                    |
  | From tag v6-08-06, 2 March 2017                            |
  | Try '.help', '.demo', '.license', '.credits', '.quit'/'.q' |
   ------------------------------------------------------------

Alex Style being used.
OBJ: TStyle     alexStyle       Alex's Style : 0 at: 0x201fcd0
root [0]
Attaching file dune1x2x6_optical_tutorial_sim_hist.root as _file0...
(TFile *) 0x26ee220
root [1] .ls
TFile**         dune1x2x6_optical_tutorial_sim_hist.root
 TFile*         dune1x2x6_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()
(Bool_t) true
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
root [4]

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:

dune1x2x6_optical_tutorial_xscan.fcl
dune1x2x6_optical_tutorial_yscan.fcl
dune1x2x6_optical_tutorial_zscan.fcl

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 dune1x2x6_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:

dune1x2x6_optical_example_yscan_gen.root
dune1x2x6_optical_example_yscan_hist.root

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, passing as an argument to Tutorial2() the number of the photon detector you want to look at. They number 0 to 119:

% root
   ------------------------------------------------------------
  | Welcome to ROOT 6.08/06                http://root.cern.ch |
  |                               (c) 1995-2016, The ROOT Team |
  | Built for macosx64                                         |
  | From tag v6-08-06, 2 March 2017                            |
  | Try '.help', '.demo', '.license', '.credits', '.quit'/'.q' |
   ------------------------------------------------------------

Alex Style being used.
OBJ: TStyle    alexStyle    Alex's Style : 0 at: 0x7f985c2bc600
root [0] .L Tutorial2.C
root [1] Tutorial2(20)

You should get a plot that looks like this:

Try a few channels and see if you can figure out where they're located based on how they respond to muons which are simulated in different locations.

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

c1->Print("yscan_pd20.png")
Now, try the following:
  • Once you've located a few of the PDs in y, try doing the scans in x and z. Again, try to locate the PDs in each dimension. (See p. 3 of this presentation if you’re having trouble.)

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 dune_pardata package. You likely don't have it checked out so to look at it directly you'll need to go to the base product:

$DUNE_PARDATA_DIR/PhotonPropagation/LibraryData/lib_35ton_v5_20150721.root

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 dune1x2x6_optical_tutorial_libanalysis.fcl. Run it and let's look at the output.

lar -c dune1x2x6_optical_tutorial_libanalysis.fcl

This produces a file called PhotonLibraryFile_dune1x2x6example_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_dune1x2x6example_hists.root as _file0...
root [1] .ls
TFile**        PhotonLibraryFile_dune1x2x6example_hists.root
 TFile*        PhotonLibraryFile_dune1x2x6example_hists.root
  KEY: TDirectoryFile    libana;1    libana (PhotonLibraryAnalyzer) folder
root [2] libana->cd()
(Bool_t)1
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 dune1x2x6_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::dune1x2x6_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 dune1x2x6_optical_tutorial_digi.fcl -s dune1x2x6_optical_tutorial_sim_gen.root

As usual, it will produce two files, dune1x2x6_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 dune1x2x6_optical_tutorial_digi_hist.root as _file0...
root [1] .ls
TFile**         dune1x2x6_optical_tutorial_digi_hist.root
 TFile*         dune1x2x6_optical_tutorial_digi_hist.root
  KEY: TDirectoryFile   daq;1   daq (SimWiredune1x2x6) 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()
(Bool_t)1
root [3] .ls
TDirectoryFile*         opdigiana       opdigiana (OpDetDigiAnaDUNE) folder
 KEY: TH1D      event_1_opchannel_0_waveform_0;1
 KEY: TH1D      event_1_opchannel_1_waveform_0;1
 KEY: TH1D      event_1_opchannel_2_waveform_0;1
 KEY: TH1D      event_1_opchannel_3_waveform_0;1
 KEY: TH1D      event_1_opchannel_4_waveform_0;1
 KEY: TH1D      event_1_opchannel_4_waveform_1;1
 KEY: TH1D      event_1_opchannel_5_waveform_0;1
 KEY: TH1D      event_1_opchannel_5_waveform_1;1
 KEY: TH1D      event_1_opchannel_6_waveform_0;1
...
 KEY: TH1D      event_1_opchannel_478_waveform_0;1
 KEY: TH1D      event_1_opchannel_479_waveform_0;1
 KEY: TH1D      event_1_opchannel_479_waveform_1;1

Each histogram corresponds to one digitized waveform. The naming convention is:

  • Event: The simulated event number (in this case, there is only 1)
  • OpChannel: The channel the waveform was read out on. In the 1x2x6 geometry there are 12 APAs, each APA has 10 photon detectors, and each photon detector has 4 electronics channels, for a total of 480 channels.
  • Waveform: It is possible for a single channel to have several digitized waveforms during an event, and this increments for each one seen.

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 dune1x2x6_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 OpHit and OpFlash modules which build the OpHits and OpFlashes:

      # photon detector reconstruction
      ophit:     @local::dunefd_ophit
      opflash:   @local::dunefd_opflash

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

      opflashana:  @local::dunefd_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 dune1x2x6_optical_tutorial_reco.fcl -s dune1x2x6_optical_tutorial_digi_gen.root

As usual, it will produce two files, dune1x2x6_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 dune1x2x6_optical_tutorial_reco_hist.root as _file0...
(TFile *) 0x26b05e0
root [1] .ls
TFile**         dune1x2x6_optical_tutorial_reco_hist.root
 TFile*         dune1x2x6_optical_tutorial_reco_hist.root
  KEY: TDirectoryFile   opflashana;1    opflashana (OpFlashAna) folder
root [2] opflashana->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
root [3]

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

PerOpHitTree

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

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: Looking at photon detector and other information using the AnalysisTree

The above trees give a lot of detail about the photon detector data, but they do not allow you to connect what is happening with the photon detectors to MC truth or other activity in the detector. To do that, we will use a different module to produce an AnalysisTree. An example fhicl file for producing analysis trees can be found in the tutorial directory:

dune1x2x6_optical_tutorial_stdana.fcl

However, it will not work when run on the previous tutorial steps since they contain only photon detector information and are missing the full, detailed reconstruction information that the tree is looking for. In this tree, each entry corresponds to 1 chunk of drift time containing a single Supernova neutrino interaction and many Ar39 interactions.

So, a file has been prepared in advance for you to use:

/dune/data/users/ahimmel/PDTutorial/dune1x2x6_optical_tutorial_ana_hist.root

This is just a simple root tree, so you can just open it up in root:

root [0]
Attaching file /dune/data/users/ahimmel/PDTutorial/dune1x2x6_optical_tutorial_ana_hist.root as _file0...
(TFile *) 0x349a1f0
root [1] .ls
TFile**         /dune/data/users/ahimmel/PDTutorial/dune1x2x6_optical_tutorial_ana_hist.root
 TFile*         /dune/data/users/ahimmel/PDTutorial/dune1x2x6_optical_tutorial_ana_hist.root
  KEY: TDirectoryFile   analysistree;1  analysistree (dune/AnaTree/AnalysisTree) folder
root [2] analysistree->cd()
(bool) true
root [3] .ls
TDirectoryFile*         analysistree    analysistree (dune/AnaTree/AnalysisTree) folder
 KEY: TTree     anatree;41      analysis tree
 KEY: TTree     anatree;40      analysis tree
 KEY: TTree     pottree;1       pot tree
root [4]  anatree->Show(0)

That will show a lot of information you need to scroll through. For the photon detectors, it only contains the information about Flashes. The flash properties are stored in vectors where the same index represents the same flash in each vector.

 no_flashes      = 444
 flash_time      = 0.411776,
                  1730.53, 1267.42, -285.403, 1742.34, 1923.21,
                  -28.9238, -20.2469, -1908.98, 940.707, -2118.62,
                  870.112, 842.709, 180.073, 1967.61, -2060.83,
                  400.432, -2006.44, -1058.97, -1894.51
 flash_pe        = 53.8876,
                  15.2468, 14.7941, 14.7006, 13.5371, 12.868,
                  12.3219, 11.2882, 11.0486, 10.9877, 10.9648,
                  10.9304, 10.8497, 10.1117, 10.038, 9.29115,
                  9.18155, 9.10811, 9.05225, 9.03688
 flash_ycenter   = 287.234,
                  -376.052, 100.919, 145.873, 51.7566, -100.555,
                  83.6191, 310.444, 243.51, -405.319, 236.438,
                  164.306, -174.435, 205.75, 420.142, -130.796,
                  -45.3386, -85.197, 333.44, 214.669
 flash_zcenter   = 690.682,
                  699.819, 931.916, 752.212, 1180.99, 721.625,
                  620.35, 197.824, 809.29, 1277.27, 732.044,
                  746.867, 734.701, 509.97, 461.001, 712.308,
                  535.99, 412.075, 763.636, 719.923
 flash_ywidth    = 161.972,
                  203.758, 277.77, 207.892, 46.2917, 227.3,
                  237.422, 163.701, 306.305, 0, 229.866,
                  266.116, 310.677, 254.111, 206.524, 278.828,
                  153.888, 126.288, 310.56, 206.799
 flash_zwidth    = 200.088,
                  313.498, 377.492, 344.503, 221.672, 251.683,
                  226.022, 177.594, 200.578, 0, 408.611,
                  266.266, 396.094, 295.252, 439.767, 422.69,
                  303.382, 456.234, 284.948, 346.569
 flash_timewidth = 0.417969,
                  0.410156, 0.476562, 0.445312, 0.472656, 0.285156,
                  0.332031, 0.234375, 0.425781, 0.269531, 0.375,
                  0.464844, 0.394531, 0.328125, 0.214844, 0.410156,
                  0.222656, 0.09375, 0.382812, 0.34375

This event has > 400 flashes, and Show only shows you the first handful of them. The flashes are sorted by PE so the first entry will often (though not always) be the flash associated with the SN event. You can draw the PE distribution of just those flashes as follows:

root [5] anatree->Draw("flash_pe[0]>>h0(200,0,200)")
root [6] c1->Print("flash_pe_0.png")
Info in <TCanvas::Print>: png file flash_pe_0.png has been created

Tutorial 7: How to match flashes to truth information

The above show's you how to get both flash and track reconstructed information in the same tree, but it is very unclear which flash corresponds to which bit of reconstructed track. This tutorial uses an example algorithm which will give you flashes and truth information together in the same tree. It also shows how to use the photon backtracker to ask about which MC truth particles contributed to a particular flash.

The fhicl to run is dune1x2x6_optical_tutorial_flashmatch.fcl and is available in the tutorial directory. It runs a module called FlashMatchAna_module.cc which is found in the OpticalDetector subdirectory of DUNE. This is a highly recommended launching off point for writing your own analyzers of photon detector information.

Running the analyzer

So, let's run this and see what we find. We will be running over some pre-made files with a mix of supernova and radiological background events in them. Standard dune datasets can be found at http://dune-data.fnal.gov. For this example we will use a file from the MCC10 production dataset called snb_timedep_radio_dune10kt_1x2x6_snb_timedep_bkg_reco. The example below shows how to retrieve the name of one file from that dataset, and then stream it from its location on PNFS. (The file is quite large, 4.9 GB, so it is better not to have to copy it to your disk space if you can avoid it).

$ setup_fnal_security
$ samweb list-files "defname: snb_timedep_radio_dune10kt_1x2x6_snb_timedep_bkg_reco with limit 1" 
snb_timedep_radio_dune10kt_1x2x6_9623_20171224T033911_reco0.root
$ samweb2xrootd snb_timedep_radio_dune10kt_1x2x6_9623_20171224T033911_reco0.root
xroot://fndca1.fnal.gov:1094/pnfs/fnal.gov/usr/dune/tape_backed/dunepro/mc/dune/full-reconstructed/02/00/13/50/snb_timedep_radio_dune10kt_1x2x6_9623_20171224T033911_reco0.root

The script samweb2xrootd creates a URL that root knows how to open to access the file. xrootd access requires some special security setup (a certificate and a grid proxy) which are set up by setup_fnal_security. You should only need to run it about once per week, and getting a Error in <TNetXNGFile::Open>: [FATAL] Auth failed error is a sign you need to run it again with the --force option.

This allows you to just pass the filename, along with that script, when running over a file with `lar`. The example below takes about 2 minutes to run:

$ lar -c dune1x2x6_optical_tutorial_flashmatch.fcl `samweb2xrootd snb_timedep_radio_dune10kt_1x2x6_9623_20171224T033911_reco0.root`
%MSG-i MF_INIT_OK:  Early 03-Aug-2018 09:42:25 CDT JobSetup
Messagelogger initialization complete.
%MSG
Info in <TGeoManager::Import>: Reading geometry from file: /cvmfs/dune.opensciencegrid.org/products/dune/dunetpc/v07_00_00/gdml/dune10kt_v2_1x2x6.gdml
Info in <TGeoManager::TGeoManager>: Geometry GDMLImport, Geometry imported from GDML created
.
.
.
TrigReport ---------- Event  Summary ------------
TrigReport Events total = 100 passed = 100 failed = 0

TrigReport ------ Modules in End-Path: end_path ------------
TrigReport  Trig Bit#        Run    Success      Error Name
TrigReport     0    0        100        100          0 flashmatch

TimeReport ---------- Time  Summary ---[sec]----
TimeReport CPU = 26.221013 Real = 100.109952

MemReport  ---------- Memory  Summary ---[base-10 MB]----
MemReport  VmPeak = 1468.03 VmHWM = 476.832

Art has completed and will exit with status 0.

Looking at the output tree

The above example produces a file called dune1x2x6_optical_tutorial_flashmatch_hist.root.

$ root dune1x2x6_optical_tutorial_flashmatch_hist.root
root [0]
Attaching file dune1x2x6_optical_tutorial_flashmatch_hist.root as _file0...
(TFile *) 0x23ab6f0
root [1] .ls
TFile**         dune1x2x6_optical_tutorial_flashmatch_hist.root
 TFile*         dune1x2x6_optical_tutorial_flashmatch_hist.root
  KEY: TDirectoryFile   flashmatch;1    flashmatch (FlashMatchAna) folder
root [2] flashmatch->cd()
(bool) true
root [3] .ls
TDirectoryFile*         flashmatch      flashmatch (FlashMatchAna) folder
 KEY: TTree     FlashMatchTree;1        FlashMatchTree
 KEY: TTree     LargestFlashTree;1      LargestFlashTree
 KEY: TTree     SelectedFlashTree;1     SelectedFlashTree
 KEY: TEfficiency       recoEfficiencyVsE;1
 KEY: TEfficiency       recoEfficiencyVsX;1
 KEY: TEfficiency       largestEfficiencyVsE;1
 KEY: TEfficiency       largestEfficiencyVsX;1
 KEY: TEfficiency       selectedEfficiencyVsE;1
 KEY: TEfficiency       selectedEfficiencyVsX;1
root [5]

The output root file contains 3 trees. All three trees have one entry per event and contain the same truth information about the original event, but contain different information about the flashes in the event.

  • FlashMatchTree -- has vectors with the properties of every flash within one drift time of the MC truth event.
  • LargestFlashTree -- has information only about the largest flash within one drift time.
  • SelectedFlashTree -- has information only about the largest flash within the distance cut defined in the fhicl (240 cm by default).

For example, this is the first entry in the SelectedFlashTree:

root [4] SelectedFlashTree->Show(0)
======> EVENT:0
 EventID         = 822401
 TrueX           = -270.76
 TrueY           = -199.284
 TrueZ           = 125.702
 TrueT           = 0
 DetectedT       = 1661.7
 TrueE           = 0.0120685
 TruePDG         = 12
 NFlashes        = 71
 FlashID         = 14
 YCenter         = -155.892
 ZCenter         = 347.709
 YWidth          = 0
 ZWidth          = 0
 Time            = 499.868
 TimeWidth       = 0
 TimeDiff        = 1161.83
 TotalPE         = 5.32417
 NOpDets         = 120
 NHitOpDets      = 1
 PEsPerOpDetVector = (vector<float>*)0x50a0420
 Purity          = 0
 Distance        = 226.208

It contains some truth information about the signal in the event as well as the properties of the largest flash that was found. Note that the Purity in this particular event = 0, so this is an event where the wrong flash was identified.

By comparison, here is the first entry in the FlashMatchTree:

root [5] FlashMatchTree->Show(0)
======> EVENT:0
 EventID         = 822401
 TrueX           = -270.76
 TrueY           = -199.284
 TrueZ           = 125.702
 TrueT           = 0
 DetectedT       = 1661.7
 TrueE           = 0.0120685
 TruePDG         = 12
 NFlashes        = 71
 FlashIDVector   = (vector<int>*)0x6b668a0
 YCenterVector   = (vector<float>*)0x6b6b060
 ZCenterVector   = (vector<float>*)0x6adbe80
 YWidthVector    = (vector<float>*)0x6ae1390
 ZWidthVector    = (vector<float>*)0x6ae8850
 TimeVector      = (vector<float>*)0x2376660
 TimeWidthVector = (vector<float>*)0x6be79c0
 TimeDiffVector  = (vector<float>*)0x6b0cbe0
 TotalPEVector   = (vector<float>*)0x6ace240
 NOpDets         = 120
 NHitOpDetVector = (vector<int>*)0x6af1600
 Purity          = (vector<float>*)0x6ade660
 Distance        = (vector<float>*)0x6ae59b0

It contains some truth information about the signal in the event, and vectors where the same entry in each vector corresponds to the same flash. Note the Signal and Purity vectors. The flash with the highest purity in the event is marked with the boolean Signal in order to make it easier to plot signal and background. For example, we can plot the PE for flashes marked as signal and compare to the others.

root [5] FlashMatchTree->Draw("TotalPEVector>>hSignal", "Signal")
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
(long long) 78
root [6] FlashMatchTree->Draw("TotalPEVector>>hBackground", "!Signal", "histsame")
(long long) 2180
root [7] hBackground->SetLineColor(kRed)
root [8] hBackground->Scale( 10./ hBackground->GetMaximum() )

Try looking at other distributions to see if they show a difference between signal and non-signal flashes.

The analyzer also produces some TEfficiency objects. Think of these objects as containers for two histograms: the numerator (selected events) and the denominator (all events). When you call Draw on a TEfficiency plot it does the division between those two histograms to show the efficiency. One key advantage of using TEfficiency objects is that if you combine root files with hadd you will still get the correct efficiency, which is not true if you had stored the efficiency in a simple histogram.

In this case, the efficiencies are plotted against two different true variables: the true neutrino energy and the true vertex position in X for 3 variants:

  • recoEfficiencyVsE/X -- Did a flash get reconstructed at all that contained some signal photons?
  • largestEfficiencyVsE/X -- By choosing the largest flash, how often did we get choose a flash associated with the signal event?
  • selectedEfficiencyVsE/X -- By choosing the largest flash within a distance cut, how often did we get choose a flash associated with the signal event?

Tutorial 8: How to resimulate MC to change photon detector efficiency

At this stage the final performance of the photon detector system isn't known. So, we need to run simulations to examine physics performance at a variety of detector performances. Recent MCC 10 simulation files were produced in such a way that we can re-run only the final "digitizer" part of the simulation and reconstruction when changing detector efficiency (up to a point), and we do not need to go back to Geant4.

You can see what efficiency was used in existing files using a tool called config_dumper. It will give you the full fhicl configuration used to produce a file -- modules run as-is and services when run with the -S option. The detector efficiency is set with the QuantumEfficiency parameter in the DUNEOpDetResponse, which is a particular implementation of the OpDetResponseServiceInterface. For DUNE simulation this is not really quantum efficiency, so much as the total efficiency of the collection system (except the parts which are dependent on position on the bar, which are part of the underlying photon library). It also includes a factor of 0.7 to account of the absorption of light on the wires and meshes which sit in front of the photon detectors but which are not included in the optical simulation. config_dumper uses a lot of output so you should use less to navigate it more easily.

$ config_dumper -S snb_timedep_radio_dune10kt_1x2x6_9623_20171224T033911_reco0.root | less
=============================================
Processing file: snb_timedep_radio_dune10kt_1x2x6_9623_20171224T033911_reco0.root
# Read SQLiteDB from file, total size: 184.0 KiB.

BackTrackerService: {
   BackTracker: {
      G4ModuleLabel: "largeant" 
      MinimumHitEnergyFraction: 1e-1
   }
}

DUNEGeometryHelper: {
}

DUNEOpDetResponse: {
   ChannelConversion: "fast" 
   FracLong: 7.1e-1
   FracShort: 2.9e-1
   LambdaLong: 225
   LambdaShort: 4.3
   LightGuideAttenuation: true
   LongAxis: "z" 
   QuantumEfficiency: 2.87e-3
   WavelengthCutHigh: 10000
   WavelengthCutLow: 0
}

You can see above that the efficiency used in 0.00287. In this tutorial we will try doubling the efficiency to 0.00574 and see how it affects the physics performance.

The fhicl dune1x2x6_optical_tutorial_resimulate.fcl makes this change (note that we make the change in the Interface, not the implementation):

########################################
# Set a new photon detector efficiency #
########################################

#                                                            Effecitve Area  Comment
#services.OpDetResponseInterface.QuantumEfficiency: 0.00287   # 4.05 cm2     This is the nominal
#services.OpDetResponseInterface.QuantumEfficiency: 0.00361   # 5.10 cm2     Improved dip-coated desing
#services.OpDetResponseInterface.QuantumEfficiency: 0.005265  # 7.44 cm2     Current double-shift
services.OpDetResponseInterface.QuantumEfficiency:  0.009059  # 12.8 cm2     Highest ARAPUCA at Tallbo
#services.OpDetResponseInterface.QuantumEfficiency: 0.01063   # 15.0 cm2     Potential improved double-shift
#services.OpDetResponseInterface.QuantumEfficiency: 0.016277  # 23.0 cm2     1.3% efficiency ARAPUCA
and then runs many of the tutorial stages we've run separately before all together in a single fhicl:
   # Run both detector simulation and reconstruction
   producers:
   {
      opdigi:    @local::dunefd_opdigi_threegang    # simple digitizer with no noise and high saturation
      ophit:     @local::dunefd_ophit
      opflash:   @local::dunefd_opflash
      rns:       { module_type: "RandomNumberSaver" }
   }

   # Run both analyzers we've looked at
   analyzers:
   {
      opflashana:  @local::dunefd_opflashana
      flashmatch:  @local::standard_flashmatchana
   }

Before running re-simulation, however, there are some checks that need to be performed related to a feature of the simulation called "scintillation prescaling." While Argon in principle produces 10's of thousands of photons per MeV deposited, we do not actually simulate all of those photons for every event. We pre-scale down the number we produce since we know the efficiency of the detectors will throw away the vast majority of those photons anyway. But, what this means is that we need to take great care with the value set for this prescaling. It needs to be the same in the input file and in the fhicl being run on that file, and it must be greater than the quantum efficiency being simulated.

A short script does the necessary checks for you is included in the Tutorial directory:

./CheckPhotonPrescale.sh <fhicl file> <input file> (or `samweb2xrootd <input filename from sam>`)

If it returns Good to go! then you can proceed with the simulation. If not, it will tell you what changes are necessary in your fhicl file in order to proceed.

After running, you will get a new histogram file, dune1x2x6_optical_tutorial_resimulate_hist.root whose contents matches those of the previous tutorials. Use those same tools to see how the signal and background have changed now that the photon detectors are twice as efficient.

To run over multiple files, you can use project.py and job submission to accomplish this. For a more complete overview on how project.py works, you can read through the tutorial here: https://cdcvs.fnal.gov/redmine/projects/dunetpc/wiki/Using_project_python

Below is an example XML script, job_dune1x2x6_optical_tutorial_resimulate.xml, which runs dune1x2x6_optical_tutorial_resimulate.fcl over many files contained in a SAM definition (you can read more about SAM here: https://wiki.dunescience.org/wiki/SAM):

<!-- XML script to run dune1x2x6_optical_tutorial_resimulate.fcl over SAM definition -->
<?xml version="1.0"?>

<!DOCTYPE project [
<!ENTITY release "v07_06_00"> <!-- Change this line to a different LArSoft version -->
<!ENTITY file_type "mc">
<!ENTITY run_type "physics">
<!ENTITY name "dune1x2x6_optical_tutorial_resimulate">
<!ENTITY tag "devel">
<!ENTITY userid "ahimmel"> <!-- Change this to an appropriate userid -->
]>

<project name="&name;">

 <!-- Group -->
  <group>dune</group>
  <numevents>10000</numevents>
  <os>SL6</os>

<resource>DEDICATED,OPPORTUNISTIC</resource>
  <memory>2000</memory> <!-- Change this line to request more memory -->

  <!-- Larsoft information -->
  <larsoft>
      <tag>&release;</tag>
      <qual>e17:prof</qual>
      <!-- Use the <local> line to specify a specific LArSoft tarball -->
  </larsoft>

  <stage name="resimulate">
      <fcl>dune1x2x6_optical_tutorial_resimulate.fcl</fcl> <!-- Change this line to a specific fhicl file, if applicable -->
      <inputdef>mcc11_snb_monoenergetic_marley_clean</inputdef> <!-- Change this line to consider a different sample; use the SAM definition! -->
      <outdir>/pnfs/dune/scratch/users/&userid;/&release;/&name;</outdir> <!-- Change this if you want to store the output somewhere outside of scratch -->
      <workdir>/pnfs/dune/scratch/users/&userid;/&release;/&name;</workdir>
      <numjobs>1012</numjobs> <!-- Change this line to split the stage into a different number of jobs -->
      <maxfilesperjob>1</maxfilesperjob>
      <datatier>detector-simulated</datatier>
      <jobsub>--expected-lifetime=3h</jobsub> <!-- Change this line to request more time -->
  </stage>

<!-- file type -->
  <filetype>&file_type;</filetype>

  <!-- run type -->
  <runtype>&run_type;</runtype>

</project>

You use this XML script to submit these jobs to the grid with project.py with the following command:

project.py --xml job_dune1x2x6_optical_tutorial_resimulate.xml --stage resimulate --submit --clean

where

--xml denotes the XML script you are using;
--stage denotes the stage within the XML script you would like to submit;
--submit actually submits the jobs to the grid; and 
--clean verifies that the output directory is empty before you submit the jobs.

Tutorial 9: How to examine and reformat raw 35ton photon detector data

First, you need to obtain some raw data to work with. Files are accessed from the SAM database, so you need to be on the Fermilab computers. The general command is:

ifdh_fetch -e dune <filename>

You may need to run setup fife_utils in order to use ifdh_fetch.

Some options you can use for files:
  • v35t_r0000101_s0001_test_raw.root
    • Vertical slice (warm), reading out LED flashes. This is what the tutorial uses.
  • lbne_r005130_sr01_20151203T025933.root
    • 35ton data with photon detector readouts taken only when externally triggered by the muon paddles (run rces_and_ssps_and_ptb).
    • Comes from this SAM query: samweb list-files "data_tier raw and lbne_data.run_mode rces_and_ssps_and_ptb and lbne_data.detector_type like %ssp%"
  • lbne_r014906_sr01_20160306T230611.root
    • 35ton data with photon detector readouts taken only also self-triggering, using the weeder to store only headers for less interesting data (run rces_and_stssps_and_ptb).
    • Comes from this SAM query: samweb list-files "data_tier raw and lbne_data.run_mode rces_and_stssps_and_ptb and lbne_data.detector_type like %ssp%"

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:

ifdh_fetch -e dune v35t_r0000101_s0001_test_raw.root
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

PROCESS NAME | MODULE LABEL | PRODUCT INSTANCE NAME | DATA PRODUCT TYPE............ | PRODUCT FRIENDLY TYPE | SIZE
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:

raw::OpDetPulses_ssptooffline_offlinePhoton_RunSSPToOffline.

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:

If looking at files that included self-triggering (the stssp configuration), then there will be an additional photon detector data product built from the header information of photon detector data where the waveform wasn't stored (generally photon detector signals which occurred when there was no TPC data). These objects are identical in structure to the OpHit object described early in the tutorial.

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.

How to make a photon library

Bonus Tutorial: How to run full photon simulation

This tutorial will not be necessary for everyone, so you probably don't need to go through it until you need it.

How to run full photon simulation

Bonus Tutorial: How to generate your own correction curves for the semi-Analytic approach

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 correction curves will only need to be remade occasionally. Also note that the tutorial has been written in a different project (SBN), but everything explained there is perfectly applicable to DUNE.

How to make correction curves for Semi-Analytic simulation mode