- Table of contents
- Photon Simulation Tutorials
- Tutorial 0: Setup your environment
- Tutorial 1: Generate photon signals
- Tutorial 2: Compare photon signals
- Tutorial 3: How to examine a photon library directly
- Tutorial 4: How to generate digitized photon detector signals
- Tutorial 5: How to run photon detector reconstruction and examine the output
- Tutorial 6: Looking at photon detector and other information using the AnalysisTree
- Tutorial 7: How to match flashes to truth information
- Tutorial 8: How to resimulate MC to change photon detector efficiency
- Tutorial 9: How to examine and reformat raw 35ton photon detector data
- Bonus Tutorial: How to generate your own photon library
- Bonus Tutorial: How to run full photon simulation
- Bonus Tutorial: How to generate your own correction curves for the semi-Analytic approach
- Bonus Tutorial: How to make computable graph for photon fast simulation
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 = 0The 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. (While this was true for the example files here, modern simulation includes the effect of the wires as part of the parameterized 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 ARAPUCAand then runs many of the tutorial stages we've run separately before all together in a two fhicls (the producers and analyzers need to be run separately for technical reasons):
Resimulate:
# 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" }
}
Resimana:
# 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 both, you will get a new histogram file, dune1x2x6_optical_tutorial_resimana_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). Note, the ScintPreScale value for MCC11 is not the same as it was for MCC10 -- make sure to check it!
<!-- 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.
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.
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
Bonus Tutorial: How to make computable graph for photon fast simulation¶
This tutorial will not be necessary for everyone, so you probably don't need to go through it until you need it.