Prog. Theor. Exp. Phys. 2013, 073C01 (20 pages)
DOI: 10.1093/ptep/ptt046

Performance of a remotely located muonradiography system to identify the inner
structure of a nuclear plant
H. Fujii1 , K. Hara2,∗ , S. Hashimoto2† , F. Ito2 , H. Kakuno3 , S.H. Kim2 , M. Kochiyama1 , K.
Nagamine1 , A. Suzuki1 , Y. Takada2 , Y. Takahashi2‡ , F. Takasaki1 , and S. Yamashita41

KEK, 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan
University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8571, Japan
3
Tokyo Metropolitan University, 1-1 Minami-Osawa, Hachioji, Tokyo 192-0397, Japan
4
University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan∗
E-mail: hara@hep.px.tsukuba.ac.jp2...............................................................................
The performance of a muon radiography system designed to image the inner structure of a
nuclear plant located at a distance of 64 m is evaluated. We conclude that there is an absenceof fuel in the pressure vessel during the measurement period and profile the fuel material placedin the storage pool. The obtained data also demonstrate the sensitivity of the system to waterlevel changes in the reactor well and the dryer-separator pool. It is expected that the system could
easily reconstruct a 2 m cubic fuel object. By operating multiple systems, typically four identicalsystems, viewing the reactor from different directions simultaneously, detection of a 1 m cubic
object should also be achievable within a period of a few months.

...............................................................................
Subject Index

1.

C50, F00, H21

Introduction

Soon after the Fukushima Daiichi reactor accident, one of the authors (K.N.) proposed that it could bepossible to inspect damaged reactors using cosmic rays. Some pioneering works [1–5] have alreadyused cosmic muons to study the inner structure of big objects such as volcanoes and hot melting
furnaces: this method is often referred to as muon imaging or muon radiography. By accumulating
millions of muon tracks traversing through such heavy and big objects, one can image their inner
structure, much like medical X-ray imaging.
A team was formed by the people of KEK, the University of Tsukuba, and the University of Tokyo
to consider the possible applications of muon imaging technology to troubled nuclear reactors. It
immediately became clear that the situation surrounding the reactors, e.g. the radiation level, would
not allow an individual to bring the detector system inside the reactor building. The radiation levels
just outside of the Fukushima Daiichi building were so high that even some protection of the detector
against radiation needed to be carefully considered. The accuracy of the track position measurement
was re-evaluated so that the detector system could provide a clear image of the inner structure even

† Now at NEDO, 1310 Omiya, Saiwai-ku, Kawasaki, Kanagawa 212-8554, Japan‡ Now at JAXA, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 252-5210, Japan© The Author(s) 2013. Published by Oxford University Press on behalf of the Physical Society of Japan.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0/),
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Received February 2, 2013; Accepted April 27, 2013; Published July 1, 2013

PTEP 2013, 073C01

H. Fujii et al.

if placed outside of the reactor building. An attempt to image the fuel in a reactor using emulsion
plates located directly underneath the reactor has been reported [6]. Setting a detector at the sameplace would be difficult in view of the accessibility required for data collection. In light of these
conditions, the KEK–Tsukuba team built a detector system using finely segmented hodoscopes consisting of 1 cm (wide) by 100 cm (long) scintillator bars. Wavelength shifting fibers were embeddedin the bars and coupled with MPPCs (multi-pixel photon counters or silicon PMTs) for readout.
The signal readout and trigger processing was performed using an FPGA (field programmable gatearray)-based system. These detector technologies, each of which had been originally developed forvarious individual physics experiments, were integrated into a muon radiography system that could
be used to evaluate the inner structure of the huge nuclear plant, a task which would otherwise bevery difficult.
In order to evaluate the performance of the detector, we carried out a test measurement by imaging
a GM MK-II type nuclear reactor at the Japan Atomic Power Company (JAPC). The system wasinstalled in March 2012 about 64 m away from the reactor center, and this report describes the designof the detector system, operational experience, and outcome of the test measurements carried out at
the JAPC.

2. Detector system
2.1. Overview of the detector
The detector system illustrated in Fig. 1 consists of four X–Y units with a 61 cm long 1.4 T permanent
magnet located between the second and third X–Y units. Each X–Y plane was constructed from two
detector planes arranged at 90◦ , each plane consisting of 100 scintillating strips of width 1 cm. The
distance between Unit-1 and Unit-2 (composing the front tracker) was 1.5 m, whilst the distance
between Unit-3 and Unit-4 was 0.75 m (compromising between the acceptance of cosmic muons andthe momentum resolution). The elevations of the four X–Y units were adjustable with respect to thetoroidal magnet.
Each X–Y unit (see Fig. 2) was equipped with a data-acquisition (DAQ) box that collected 200
MPPC signals; 100 each of X and Y. The X–Y coincidences (after applying thresholds to the signals
in a given time window) were transferred via a LAN to a PC, where the trigger decision was made.
2/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 1. The detector, consisting of four 1 × 1 m X–Y (scintillator bar) units and a toroidal magnet. The DAQ
box looks for a coincidence in the X–Y unit and sends the data to the PC. The clocks of the four DAQ boxes
were synchronized by a clock generator module.
These detector technologies, each of which had been originally developed for various individual physics
experiments, were integrated into a muon radiography system that could be used to evaluate the inner structureof the huge nuclear plant, a task which would otherwise be very difficult.

PTEP 2013, 073C01

H. Fujii et al.

The timing of the four DAQ boxes was reset every 1 µs using a common clock generator that provided
NIM-logic signals to all four of the DAQ boxes. The whole system was housed in a 20 ft shipping
container, reinforced with steel bars and equipped with air conditioning, which weighed 3.5 t. The
total weight of the system amounted to 10 t, including the toroidal magnet (5.3 t) and four X–Y units(50 kg each).

2.2.

The X–Y scintillator units

The eight detector planes, each consisting of 100 scintillating bars, were identical in design. Thebase element consisted of scintillator [7] bars of width 9.7 mm, height 6 mm (nominal), and length
1000 mm. We machined a 2 mm deep straight groove where a 1 mm φ wavelength shifting fiber(Kuraray Y11 200 ppm doped double clad plastic fiber [8]) was embedded. After cutting to a lengthof 1010 mm, both ends of the fiber were polished and one end was aluminized. The aluminum
mirror was reinforced with white paint (TiO2 -based water emulsion paint, Eau De Coat G Eco,
Nihon Paint Co., Ltd.), and the three sides of the bar were also painted white. Twenty such bars
were placed with the grooved faces set against a jig plate (where 20 ridges were machined at a
spacing of 10 mm), such that the bars were arranged at a spacing of precisely 10 mm. Five setsof such 20-bar arrangements were put in a light-tight tray, made of 8 mm thick black PMMA
bars at three sides, with a 0.5 mm thick aluminum sheet screwed to the PMMA bars. A 0.02 mm
thick white PET sheet and a 0.1 mm thick black vinyl sheet were used for light reflection andlight shielding, respectively. In order to keep the shape and flatness of the tray, four 20 × 30 mm
rectangular aluminum tubes were screwed in a direction perpendicular to the scintillator bars (seeFig. 2).
For the fiber end-plate, one hundred holes were machined at a spacing of 10 mm to a 10 mm thickblack PMMA end-plate, where wave-shifting fibers were inserted and glued. After placing the fibers
into the scintillator grooves and setting the white PET and vinyl sheets, a 2 mm thick foam sheet wasinserted, after which the tray was covered with another 0.5 mm thick aluminum sheet. The foam wasinserted for the purpose of absorbing any thickness variation of the scintillator bars. The open sideof the tray box was enclosed by screwing the fiber end-plate against a pair of enforcement bars. The
end-plates were machined with screw holes and five dowel pins, to which an MPPC array plate wasthen positioned and screwed.
3/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 2. Photograph of one X–Y unit. The DAQ box collects 200 signals from the X and Y planes. The inset isan enlarged view of the printed boards, each serving a connection of 10 MPPCs to the DAQ box via 10 coaxial
cables.

PTEP 2013, 073C01

H. Fujii et al.

Fig. 3. (Left) Distribution of bias voltages, Vb , at a gain of 7.5 × 105 , and (right) dark count rates above athreshold of 0.5 PEs. Data were taken at 25 ◦ C for all 840 delivered pieces.

The MPPC array plate

The multi-pixel photon counter (MPPC) is a product of Hamamatsu Photonics [9]; it has 667 pixels
covering a 1.3 mm square area and a typical gain of 7.5 × 105 at 71 V bias. The same devices weredeveloped for the T2K near detector [10] located at KEK. Figure 3 shows the bias voltage distribution
at a gain of 7.5 × 105 , and the dark count rates for a 0.5 photoelectron (PE) threshold for all of the 840
MPPCs used here. The devices showed quite uniform characteristics. The dark count rates decreased
to 38–86 kHz at a threshold of 1.2 PEs. The MPPCs were packaged in 6 mm φ ceramics with thephoto-coupling face covered by a 0.3 mm thick silicon resin. Each of the wave-shifting fibers wascoupled to one MPPC so that 100 MPPCs were aligned at a spacing of 10 mm along the 5 mm thickblack PMMA plate. The MPPC array plate was assembled by first placing each MPPC in drilledholes (see Fig. 4). The MPPC anode and cathode leads were then soldered onto a printed circuitboard. One circuit board, serving 10 MPPCs, was mounted with 10 2-pin headers (Hirose ElectricCo., Ltd., HNC Series); see the inset of Fig. 2. The signal was transferred via coaxial cables coupledindividually to the two leads of the MPPC, hence the signal pulse was overlaid on a typical bias
voltage of 71 V.

2.4.

The MPPC DAQ system

The MPPC data-acquisition (DAQ) system was fabricated by Bee-Beans Technology [11]. Each system dealt with 200 analog signals for the two planes of one X–Y unit, and also provided bias voltages
to the MPPCs themselves. A block diagram of this system is shown in Fig. 5. After decoupling the DC
bias voltage, the pulse signal was transferred to an LTC6409 fast differential amplifier with its gain
set to 100. The amplified signal was then discriminated using an LT1720 functioning in zero-crossing
mode, where the baseline (i.e., threshold adjustment) is provided by an individual digital-to-analog
converter (DAC). Twenty amplifier and discriminator stages were grouped, and assembled on a subboard. Ten such sub-boards were plugged into a motherboard furnished with FPGAs (for coincidence
logic generation and SiTCP communication). The motherboard also controls UDP commands to the
DACs via the LAN port. The system ran on a 125 MHz clock cycle, providing a shortest coincidence
window of 8 ns. The FPGA looks for a single cluster in both the X and Y planes in a given time
window (adjustable in units of 8 ns up to 1 µs). This cluster is defined as hits in consecutive channels
from ‘one’ up to ‘one to four’ (selectable); we chose ‘one to three’ consecutive hits in order to allow
for inclined tracks. Once this coincidence condition was met, the hit positions (in X and Y) and the
4/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

2.3.

PTEP 2013, 073C01

H. Fujii et al.

Fig. 5. Schematics of the DAQ boxes that handle 100 X and 100 Y analog signals. The bias voltage, controllable via DAC70, is common to two neighboring channels whilst the threshold is adjustable via DACofs0 and
DACofs1 individually. The DAC settings and data readouts were performed through a single LAN port.

timestamps were buffered whilst waiting to be transferred. The timestamp clock was reset and thefour DAQ boxes in use were synchronized by external pulses, provided by the clock distributor [12]that was operated at 1 kHz.
The signal line lengths in the DAQ box were designed to be identical so that the propagation inside
did not bias the coincidence judgment. The propagation time difference from the scintillator to the
5/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 4. Photo-coupling side of the MPPC array plate. The MPPC ceramics were fitted into holes of 6 mmdiameter. The innermost rectangles are MPPC sensors of 1.3 mm square.

PTEP 2013, 073C01

H. Fujii et al.

DAQ box input was 10 ns at maximum and we verified that a coincidence window of 16 ns wasappropriate for the X–Y units. The unit coincidence data (X and Y hit positions, and timestamps)
were transferred to a PC, where they were examined for coincidence in Unit-1 and Unit-2 over a 32 ns
time window. If the event was judged acceptable, these unit coincidence data were stored together
with those from Unit-3 and Unit-4 that had the closest timestamps.

2.5.

Response uniformity of the detector

2.6.

Momentum measurement

Measurement of the muon momentum provides additional information on the remote mass distribution and should help improve the radiography image by choosing tracks in an appropriate
momentum range. Also, the muon momentum information was indispensable in simulations to predict the obtained image. In light of this, we constructed a permanent toroidal magnet with a magnetic
field of 1.4 T in collaboration with NEOMAX Engineering [13]. As the results using the measured
momentum information are beyond the context of this paper and are to be reported in a separate
paper, we will not discuss this issue further here.

3. Test experiment at the JAPC
3.1. Detector installation and operation
The system was transported to the JAPC on 24 February 2012 and installed outside the reactor building at a distance to the reactor center of 64 m at an angle of 22◦ with respect to one of the building’s
walls. Although the detector started data acquisition on the same day, some detector threshold adjustments were carried out for the first few days. In the period up to 18 June 2012 (referred to here as
PERIOD-1), the detector system was aimed at the center of the fuel loading zone by adjusting therelative heights of Unit-1 and Unit-2. In the period from 28 June to 5 October (referred to here as
PERIOD-2), the elevation of Unit-1 was raised, thus aiming the detector at the fuel storage pool.

3.2.

The coincidence time distribution

The time window to judge an X–Y coincidence inside the same unit was fixed to 16 ns. The timedifference between Unit-1 and Unit-2 was then used to judge whether the event was to be recorded
or not. As seen in Fig. 6 where the time difference, T(Unit-1)−T(Unit-2), is plotted, the mean of thetime difference distribution was shifted from zero, reflecting that tracks coming from the front were
dominant in this configuration. This asymmetry was caused by the different zenith angle coveragefor the tracks coming from the front and back (Unit-1 was elevated by 20 cm in PERIOD-1), and by
the magnet behind. Events within a 32 ns time window were retained.
6/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

The light yield of the scintillator/fiber–MPPC system was typically 15 photoelectrons for β-rays
penetrating near the scintillator end opposite to the readout end, which is comfortably high enoughto detect muons. The threshold voltages of all 800 channels were first adjusted to discriminate thenoise contribution. The response uniformity in a plane was then examined using muons recorded
as coincidence in the X–Y unit placed horizontally. We accumulated about 10k events per channel
where events with only one hit in each plane were accepted. The rms variation of the counting uniformity ranged from 16% to 19%. The obtained non-uniformity includes the effects of cross talk toneighboring scintillators.

PTEP 2013, 073C01

H. Fujii et al.
3

×10
1600

hist_Dtime12

1400

Entries 6170701
Mean-2.086
RMS
11.26

1200
1000
800
600
400
200
0
-40

-30

-20

-10

0

10

Fig. 6. Time difference between Unit-1 and Unit-2 in PERIOD-1. The time was digitized in units of 8 ns.

Fig. 7. (Left) Number of recorded events versus the number of measurement days, and (right) trigger rate
calculated every 30 min. The data were accumulated over two periods, PERIOD-1 and PERIOD-2.

3.3.

Coincidence rates

Figure 7 summarizes the number of recorded events and the trigger rates during PERIOD-1 andPERIOD-2. In total, 6.3 million events and 7.1 million events were accumulated in PERIOD-1 (over114 days) and PERIOD-2 (over 99 days), respectively. In both periods the accumulated number of
events increased linearly with the number of days (except for the initial threshold adjustment period
in PERIOD-1), verifying that the detector was operated stably over the 7-month period. Each dataacquisition run lasted 30 min and the trigger rate in Fig. 7 is calculated for each run. In PERIOD-1,
the trigger rate showed a peak at 0.85 Hz, and a tail in the lower frequencies.
The MPPC gain is temperature dependent and decreases typically by 2.2% per degree. The tail wasa result of this gain variation against a fixed threshold, which happened to be unexpectedly wide due toa large variation in the room temperature from 10 ◦ C to 28 ◦ C, caused by a lack of thermal insulation
against the container walls. Towards the end of PERIOD-1, thermal insulation of the container was
accomplished. Also, at the beginning of each data-acquisition run, the room temperature was used toscale the threshold values for the MPPCs. As is evident, the trigger rate was very stable in PERIOD-2,
averaging 0.9 Hz.
The reduced trigger rate in PERIOD-1 was associated with non-uniform detector efficiency inaddition to the response non-uniformity described in Sect. 2.5. Such non-uniformity, though, has
7/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

20
30
40
T(U1) - T(U2) [ns]

PTEP 2013, 073C01

H. Fujii et al.

little impact on the image of the reactor 64 m away. This is understandable since the track projections
with the same differences of hit scintillator numbers in the two tracking units are localized within
1 m, and therefore the reactor image where we employed a bin size of 40 cm is obtained essentially byintegrating over the local non-uniformity of the detector response. The above estimate was verified
by examining the projected image using tracks coming from behind; no structure other than theacceptance was created since no massive structure was present. We also have verified that these low
trigger rate runs did not introduce any spurious projections of tracks. After these studies, we only
selected runs with a trigger rate above 0.5 Hz for further analysis.

3.4.

Track distributions

3.5.

Two-dimensional image

Two-dimensional images of detected muon tracks extrapolated to the reactor center are illustrative
to demonstrate the ability of muon radiography. As obtained images, though, they are affected by
the massive floors that prevent imaging of detailed structure of the fuel region that we are interested
in. We present a two-dimensional image, as shown in Fig. 9, where the visibility of the image was
improved by subtracting such slowly changing components. In practice, we expanded the horizontal
distributions of each vertical bin to a Fourier series and ignored the low-frequency components, the0th to 3rd series of the expansion, on reconstruction of the image. Through this procedure the effectof the geometrical acceptance is also suppressed.
In the image, the darker parts correspond to massive structures. The conical-shaped containment
vessel is visible together with the cylindrical reactor pressure vessel (RPV) at its center. Also, theright-side building wall is clearly seen, as it was almost projective to the detector.

3.6.

Sensitivity evaluation using the data

Evaluation of the detector sensitivity using only the available data is described in this section. We
present the ability to detect variation in the pool water level and the determination of the spatial
resolution of the detector system.
3.6.1. Detection of water-level variation. To investigate the possible variation of the reactor plantconfiguration, we first examined the distribution uniformity on a weekly basis, and found a distinct
variation of the event yields in specific regions. Figure 10 plots the event yield in a 4 m (horizontal)
by 9 m (vertical) region, corresponding to the dryer-separator pool, divided by that in the same 4 mhorizontal band, shown on a weekly basis. The event rate was systematically high from Week 3 to
Week 8, with transitions seen at Week 2 and Week 10. According to the work record provided after
this study was made, the water level in the reactor well had changed during this period: reducing
to −7 m on 9 March (Week 2), increasing up to −5 m on 27 April (Week 9) and −2.6 m on 1 May(Week 10), then back to the nominal 0 m on 6 May (Week 10). This variation in the water level (WL)
explains the transitions seen in Fig. 10.
8/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Since the DAQ system selected only events with a single cluster of hits per plane, tracking wasstraightforward. The reconstructed tracks obtained 64 m away from the reactor center were extrapolated to the reactor center plane (defined as Z = 0). The plots in Fig. 8 show the track distributions
along the x-axis for given slices in the y-axis at Z = 0, alongside the corresponding reactor structures
available to the public [14].

PTEP 2013, 073C01

H. Fujii et al.

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 8. Horizontal track distributions projected at the reactor center for 4 m vertical slices. The vertical regions
cover the areas (from bottom to top) just below the fuel loading zone, the fuel loading zone, just above the fuelloading zone, and further above. The recognizable dips correspond to the right-side building wall, containment
vessel walls, and pressure vessel support walls. Because of the angle of the detector with respect to the reactor
building, the projection of the left-side wall was wider in X , hence no corresponding dip was recognizable.

We were also informed by the JAPC that the well shield plug was ‘in’ between 20 March and
21 April and ‘out’ for the other periods. Whilst the plug was ‘in’, the WL reduced to −7 m.
Figure 11 shows a two-dimensional ratio of the event yield between the ‘plug-in’ period and the
9/20

PTEP 2013, 073C01

H. Fujii et al.

Fraction of events

0.22
0.21
0.20
0.19
0.18
0.17
0.16
0.15
0.14
0

2

4

6

8

10 12 14 16
Weeks since 27/2/12

Fig. 10. Event yield in the 4 × 9 m region corresponding to the dryer-separator pool, divided by the event yieldin the same 4 m horizontal band. The x-axis shows the weeks passed since 27 February 2012.

‘nominal’ period (plug was out, and WL = 0 m). The upper-right part of the figure corresponds tothe dryer-separator pool. The most distinct variation in the figure was due to a change in the waterlevel of the dryer-separator pool. The low ratio part at the top of containment vessel was due to theplug shield (1.7 m thick concrete), and the higher ratio part below the shield was due to the reduced
WL in the reactor well. Since the detector was located at an angle to the building, the same heightsof the WL in the dryer-separator pool and the reactor well did not manifest as horizontal features inthis elevation view, as can be seen from the overlaid floor structure in Fig. 11 and the density–length
distribution plot given in Sect. 4 (Fig. 15).
Figure 12 shows the vertical yield distributions for the periods with different WLs for the dryerseparator pool, normalized by the distribution obtained when WL = 0 m (5 May to 16 June). The
data points in the horizontal range corresponding to the pool region were added. The widths of thebumps reflect the reduced water level. Note that the data period was only 5 days for two of the four
periods. The dryer-separator pool was 12 m thick along the muon path, which is equivalent to a 4.8 m
10/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 9. Image obtained during PERIOD-1. The massive structures are represented by darker parts.

PTEP 2013, 073C01

H. Fujii et al.

Normalized event yield

1.5

normalized by 7/5–16/6 (WL=0 m)

1.4
1.3
1.2
1.1
1.0
0.9
27/2–8/3 (WL=0)

0.8
0.7

9/3–27/4 (-7 m)
28/4–1/5 (-5 m)
2/5–6/5 (-2.6 m)

0.6-40 -35 -30 -25 -20 -15 -10 -5 0
5Vertical [m]Fig. 12. Vertical yield distributions for periods with different water levels, normalized by the distribution at thenominal water level (WL = 0, where the vertical origin is set). The horizontal range covers the dryer-separator
pool.

thick fuel material (for a density of 2.5 g/cc). The observation of the WL changes can be interpreted
as the detection ability of the corresponding thick fuel material in the corresponding period.
3.6.2. Spatial resolution. The spatial resolution is an important parameter if we are to locate thefuel material. Here we use the data distributions to evaluate the spatial resolution. The horizontal yield
distribution across the dryer-separator pool is plotted in Fig. 13, where the distribution at WL = 0 m
was divided by that at a reduced water level (WL < 0); the dip observed here is due to the presenceof water. The dip distribution was compared to the expected distribution where the spatial resolution
11/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 11. Normalized ratio of the event yields whilst the reactor plug shield was ‘in’ and WL = −7 m (at31 days), to the event yields whilst the plug shield was ‘out’ and WL = 0 m (at 45 days). Some of the plantstructures are overlaid in black, and data points with poor statistics are excluded. A full description of the color
code is given online arXiv:1305.3423.

H. Fujii et al.
Normalized event yield

PTEP 2013, 073C01
1.6

data
expectation with σx=0.5m

1.4
1.2
1.0
0.8
0.6
0.4
-40

-20

0

20
X (m)

Event yield (No Wall-Wall)

Event yield

Wall (Data)
2800
No Wall (Data)
2600
2400
2200

700
600
500
400
300
200
100

2000

0
1800
-100
1600

-200
-4

-3

-2

-1

0

1

2

3
X [m]

-3

-2

-1

0

1

2
X [m]

Fig. 14. (Left) Horizontal event yield distribution at a vertical range including the wall under investigation (black) in comparison to that below where the wall was not included (red triangles). (Right) Flat-topped
Gaussian fit of the difference between the two curves.

was assumed to be 0.5 m, and showed good agreement, certifying the validity of this assumption forthe spatial resolution.
Another evaluation of the spatial resolution was obtained by looking at the width of a wall. One of
the walls supporting the dryer-separator pool was 1.25 m wide and was viewed as an almost projective
plane by the detector. The event yield distributions in the vertical bands both with and without thewall are shown in Fig. 14. By taking the difference between the two distributions and applying a
Gaussian fit, we obtained a spatial resolution of 56 cm with a wider wall width of 1.31 m. Since thewall was 13 m thick along the muon path and thus not completely projective, the difference obtained
from these data was then fitted to a flat-topped Gaussian, as shown in Fig. 14. This fit resulted in a
spatial resolution of 49.1 ± 2.7 cm.
12/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 13. Horizontal event yield distributions at the dryer-separator pool height, where the distribution at normal
water level (WL = 0 m) was normalized by that at a reduced water level (WL < 0). The horizontal origin isset at the center of the pool. The overlaid curve indicates the expected results where the spatial resolution was
assumed to be 0.5 m.

PTEP 2013, 073C01

H. Fujii et al.

Fig. 15. (Left) Simulated density–length distribution of the reactor. The scale is in m.w.e. (meters water equivalent). The outline of the fuel is highlighted. (Right) Evaluated density distribution inside the RPV. The density
was consistent with 1 g/cc over the entire fuel loading zone.

A simple method to evaluate the material along the muon path was developed using the density–
length method (the product of the density and the path length). The assumptions made were:(1) The observed muon flight direction was identical to the direction of the original muon entering
the structure (no scattering).
(2) The muon absorption was determined from the sum of the density–lengths of the materials
along the path. (The absorption was an exponential function of the density–length sum.)(3) The original muon intensity was dependent only on the zenith θ and azimuth φ angles. The
intensity was uniform across the azimuth angles, whereas the zenith angular dependence of
the intensity was assumed to follow cos2 θ .
(4) There was a uniformly distributed background that was independent of the zenith angle.
We incorporated these assumptions in the following equation: 
ρi ti (θ, φ)/λ + n bk
Nobs (θ, φ) = Nμ (θ, φ) · exp −where Nobs (θ, φ) is the observed muon intensity in the direction (θ, φ), Nμ is the original muon flux
entering the structure, λ is the effective absorption length, ρi is the density of the ith material, andti (θ, φ) is the length of the ith material along the muon path. In practice, the angles were defined by
the horizontal and vertical coordinate differences (l, m), and their acceptances were not uniform in(θ, φ) due to geometrical effects and muon flux. The acceptance functions f (l, m) and g(l, m) interms of the direction bin were calculated by a Monte Carlo method so that the above relation was
reduced to:Nobs (l, m) = a · f (l, m) · exp[−t (l, m)/λ] + b · g(l, m)
where Nobs is the number of muons observed in the direction bin (l, m), f (l, m) is the cos2 θ dependent acceptance function, and g(l, m) accounts for the pure geometrical acceptance of the
background.
Figure 15 shows a two-dimensional density–length distribution calculated by this method. The fuelwas assumed to be loaded in the RPV. In addition to the containment vessel and RPV, some verticalwalls and horizontal floors are recognizable.
The density distribution inside the RPV was determined by comparison with the observed event
yield, while fixing the density–lengths of the other structure components relying on the input values.
13/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

4. Comparison of the experimental results with simulations
4.1. The density–length method

PTEP 2013, 073C01

H. Fujii et al.

In this calculation, the ratio of the coefficients a and b was fixed, leaving only the normalization as afree parameter. Since the variation in absorption can be approximated as a linear function for small
density variations, the event yields for common RPV densities of 0 and 3 g/cc were calculated in a
40 × 40 cm mesh. These data-sets were used to interpolate to estimate the density in that mesh (whichis uniform along the muon path) by comparison to the observed yield. The evaluated density distribution is presented in Fig. 15, which shows that the fuel loading zone was consistent with a densityof 1 g/cc. There were some higher density parts in the upper half of the RPV that we anticipate were
due to contributions (such as scattering at the floor and the wall) that were not completely evaluated.

4.2.

Monte Carlo methodsn(P) = 3.63 · (P/(GeV/c))−0.505 .
Processing of GEANT events is time-consuming as it requires an optimization of event production.
We carried out a study to generate events randomly at Z = 0 that were then traced from a position
behind the reactor building. The scattering angle distribution for the events that were accepted by
the detector was found to be within 0.03 rad. To summarize, we first picked up the momentum of
the muon according to Jokisch et al., then generated events randomly at Z = 0 with the zenith andazimuth angles given in the way described above. This generation was repeated until the track angle
with respect to the detector center was within 0.03 rad. Successful events were traced through thebuilding and reactor material. This event production method substantially improved the reproductionof the reconstructed track distributions, compared to a method where only the tracks that would havepassed through the detector were traced (where the effects of scattering-in were not considered). As
we describe below, however, complete reproduction of the event distribution was not achieved, andwill require further studies.

4.3.

Comparison with the GEANT distribution: absence of fuelThe primary objective of the test experiment at the JAPC was to verify whether the fuel could
be imaged with our system. Throughout the experimental periods the fuel was not loaded; therefore, direct verification was not possible. Instead of relying on the GEANT simulations, the
14/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

The GEANT simulation [15] was employed to simulate the muon interactions in the reactor and thebuilding structure. Since the available plant information was limited, the building walls and floors;
the containment vessel and RPV; and the pools were the main components implemented in the simulation. The thicknesses and locations were taken from a drawing [14], where it was assumed thatthey were concrete with a density of 2.5 g/cc, except for the RPV wall and the water. The RPV wasassumed to be either empty, filled entirely with water, or filled with water and with fuel material
loaded.
The GEANT simulation traces a muon of given momentum and direction through a given material. Therefore, every critical material component needs to be taken into consideration, and reliablemuon momentum and direction distributions are indispensable. The momentum distribution was
taken from Jokisch et al. [16], who measured the distribution at a 75◦ zenith angle. For the zenith
angular dependence of the muon yield, we employed the data of the Okayama group [17], who measured the cosmic muon flux at sea level for a zenith angular range from 0 to 80◦ . They provided thezenith angular dependences for three muon momentum bins of 12.5, 22.5, and 45 GeV/c. Their zenith
angular dependences were fitted to a cosn θ function, the resultant n determined as 1.09, 0.76, and0.50, respectively. These numbers were used to obtain the momentum (P) dependence of n(P) as:H. Fujii et al.
Event yield

PTEP 2013, 073C01
8000|X|<2 m (MC SB corrected in 5<|X|<10 m)
7000

6000

5000

4000

3000
Data
2000

noWater/noFuel(MC)Water(MC)0Fuel(MC)

-8

-6

-4

-2

0

2

4

6

8
10
Vertical [m]Fig. 16. Vertical event yield distributions in the range |X | < 2 m about the fuel loading center, where GEANT
Monte Carlo (MC) distributions were corrected for sidebands (SB) of X ranging from 5 < |X | < 10 m. The
RPV was filled with water (except the noWater/noFuel case) in the vertical region from −5 m to 11 m. The
vertical origin is set at the fuel loading zone center.

images near the fuel loading zone were used to calibrate the GEANT distribution and we thenextend this understanding to the fuel region. For the GEANT simulation, three typical cases wereconsidered:(1) The pressure vessel RPV was empty (no water and no fuel). Since no material (such as fuelbrackets) was considered, this situation is equivalent to completely empty.
(2) The RPV was filled with water. Since material other than the fuel material can be approximatedby water, we expect to observe distributions obtained in this situation.
(3) A 4 m cubic fuel source was loaded. The fuel material was assumed to be UO2 , with an averagedensity of 2.5 g/cc [18]. The rest of the RPV space was filled with water.
4.3.1. Vertical distribution. The vertical event yield distribution around the fuel loading zone wasinvestigated at first. The GEANT distributions as obtained were limited in describing the detailed
structures of the plant, showing a maximum discrepancy of about 10% in the shape in the verticalrange from −8 m to +10 m about the fuel loading zone center.
The incomplete implementation of material conditions and the possibility that an inappropriate
muon momentum assumption was used were examined by comparing the experimental data andthe GEANT distributions in an area outside of the area under investigation. Figure 16 shows thevertical event yield distributions in the fuel loading range |X | < 2 m where the GEANT distributions were normalized in each vertical slice (0.4 m bin) to the corresponding data using the yields
in the horizontal sidebands of X over the range 5 < |X | < 10 m. The overall GEANT event yields
were normalized to the data in the vertical range Y < −5 m where assumptions about the waterand fuel conditions do not influence the vertical distribution. The corrected GEANT distributions
described the observed distributions fairly well, with the observed data distribution best describedby the “water-only” assumption (2).
15/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 20141000

H. Fujii et al.
Event yield

PTEP 2013, 073C01

5000
4500
4000
3500
3000

Data
noWater/noFuel(MC)Water(MC)
Fuel(MC)

2500
2000
-10

-5

0

5

10
X [m]

4.3.2. Horizontal distribution. Figure 17 shows the horizontal event yield distributions for the
4 m height range corresponding to the height of the fuel when loaded. The data distributions arecompared to the three cases employed in the GEANT simulations. A constant value (about 10% of
the event yield) was added to the simulation data to reproduce the structure in the range from thecontainment vessel wall to outside of the RPV. As shown, the characteristic distribution around thecontainment vessel was well reproduced by the simulations, which is encouraging with respect tointerpolation of the simulation distribution to inside the fuel loading zone. The data points inside
the RPV were in good agreement with the prediction where the RPV was filled with water, and can
exclude the other two cases.
Figure 18 shows similar distributions in other Y slices measured across 2 m bands. The simulated
distributions were generally in agreement with the experimental data, thereby concluding that thismethod can reliably confirm the absence of the fuel in the fuel loading zone. Note that the observed
events in [Yc + 2m, Yc + 4m] exhibit a systematical reduction around X = 2 m. This is suspected
to be caused by unaccounted-for reactor equipment located behind on the floor. The details of thematerial are not available and hence are not implemented in the simulation.

4.4.

Sensitivity evaluated using the GEANT simulations

As demonstrated in Fig. 17, the GEANT simulations can describe the observed data around the fuelregion after addition of a constant. The simulation study was extended to investigate the sensitivityof this technique for smaller fuel sizes. Figure 19 shows the horizontal distributions for different fuelsizes: 1 m to 4 m (cubic), with the distribution for the “water-only” assumption subtracted. In the
figure, we have assumed the same data-taking configuration as PERIOD-1 (114 days, at a distanceof 64 m), with the y-axis representing the number of events reduced by the existence of fuel material
in the fuel loading zone center. The vertical slice range was fixed at 1.6 m throughout.
The observation sensitivity S is expressed as:
S=[NMC (0) − NMC (w)]
√
Nobs
16/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 17. Horizontal distributions in the 4 m height range about the fuel loading center. The dips at |X | = 10 m
correspond to the containment vessel walls. The experimental data (black points) are compared with GEANT
simulations under one of three assumptions: no water and no fuel in RPV (red dotted); RPV filled entirely
with water (green solid); and with the fuel loaded with the remaining part filled with water (blue dashed). The
simulations assumed that the fuel had dimensions of 4 m × 4 m × 4 m and an average density of 2.5 g/cc.

H. Fujii et al.
Event yield / 0.4 (m)PTEP 2013, 073C01

2000[Yc+2m,Yc+4m]

1500[Yc,Yc+2m]

1500[Yc-2m,Yc]

1500[Yc-4m,Yc-2m]
Data
noWater/noFuel(MC)Water (MC)Fuel (MC)

1000
-10

0

-5

5

10

X (m)

Event Reduction/0.4 [m]Fig. 18. Horizontal distributions over a 2 m height range above and below the fuel loading zone. The height
ranges are shown with respect to the central height of the fuel loading zone, e.g., [Yc–4m,Yc–2m] refers to arange from 4 to 2 m below the fuel loading center.
200

0

-200

-400

1.0 m cubic
1.3 m cubic
2.0 m cubic-600

4.0 m cubic

-800
-10

-5

0

5

10
X [m]

Fig. 19. Horizontal yield distributions for a 1.6 m Y-slice for different fuel sizes, with the distribution obtained
for the “water-only” assumption subtracted. The fuel was assumed to be cubic with a side length varying from
1 m to 4 m and a density of 2.5 g/cc.

where Nobs is the number of observed events taken from the experimental data of PERIOD-1, and
NMC (w) is the number of simulated events for each fuel size w (0 refers to water-only). By fixing thevertical range to 1.6 m for all cases, we counted the number of events for different horizontal ranges
and calculated the sensitivity S.
We assumed that the existence of the fuel can be confirmed for S > 5, and the number of days
to achieve this condition is shown in Fig. 20 as a function of the fuel size. The quoted uncertainties are due to the variation in the side band normalization. Observations of 2 m (cubic) or larger
fuel objects should be achievable in two weeks. Also shown in the figure are the results for the
17/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

1500

H. Fujii et al.

Days to identify

PTEP 2013, 073C01
103
PERIOD-1 configuration
Fourfold statistics increase

102

10

1

0

1

2

3
4
5
Side length [m] of cubic fuelcase of a fourfold increase in event statistics that is achievable should the system be placed closer(32 m) or should four detector systems be operated simultaneously at a distance of 64 m. In thisconfiguration, a 1.3 m (1 m) fuel object should be identifiable over a one (three to four) month
period.
4.4.1. Estimation of the fuels in the storage pool. In PERIOD-2, the detector was directed towards
the fuel storage pool. The horizontal event yield distribution over a 4 m band across the pool iscompared here to the simulations. In the simulations, the fuel object was assumed to be arranged in
one rectangular block 4 m high, 10 m wide, and 0, 4, or 8 m deep. The height was a standard value
for the fuel object. The width was set wide enough so that the data distribution in the correspondingarea was not represented by the depth = 0 m distribution.
The data taken over 99 days and the simulation distributions are compared in Fig. 21. The fuelthickness as a function of the horizontal position was estimated by linear interpolation of the simulation distributions for the three thicknesses. As shown in the figure, the estimated fuel spans about
6 m horizontally with a typical thickness of 7 m (there may be a thinner part on the right). The present
estimate of the thickness relies on the chosen fuel density of 2.5 g/cc, hence the actual thickness willvary inversely with respect to the actual average density.

5.

Conclusion

We have constructed a muon radiography system consisting of 1 × 1 m tracking planes utilizing 1 cm
wide scintillator bars. The system performance was evaluated by identifying the inner structure of
the Japan Atomic Power Company (JAPC) nuclear plant. The detector system was installed 64 m
away from the center of the reactor and collected 6.3 million muons during PERIOD-1 (114 days),and 7.1 million muons during PERIOD-2 (99 days), whilst targeting the fuel loading zone and thefuel storage pool, respectively.
The detector can reconstruct an image of the characteristic reactor structures, such as the containment vessel, pressure vessel, floors, building walls, and storage pools. The image distribution wasconsistent with the case where the pressure vessel was filled with water or equivalent, and no nuclear
18/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Fig. 20. Number of days required to confirm the existence of the fuel as a function of the side length of the
cubic fuel. The data show the measurement conditions at the JAPC for 64 m from the reactor center with one
detector system (filled squares), and a fourfold increase in the statistics achieved either placing the detector
31 m from the reactor center, or operating four systems simultaneously (open squares).

H. Fujii et al.

4000

3500

3000

2500

(Data)
0m
4m
8m

2000

1500
-12

-10

-8

-6

-4

-2

0

2

4

6

Fuel thickness [m]

X [m]
10
8
6
4
2
0
-2
-4
-6

-4

-2

0

2

4

6

X [m]

Fig. 21. (Top) Estimated profile of the fuel object in the storage pool [14]. (Middle) The event yield distributions in a 4 m vertical band across the fuel storage pool (closed squares) compared to simulations where the fuelobject thickness was assumed to be 8 m (plus signs), 4 m (circles), or 0 m (stars). (Bottom) The fuel thickness
was evaluated as a function of the horizontal position (where the density was assumed to be 2.5 g/cc). The
horizontal origin is set at the center of the pool.

fuel was identified in the pressure vessel, confirming that the reactor was under maintenance during
the data-acquisition periods. We also measured the profile of the fuel material located in the fuelstorage pool. The data obtained were used to examine the sensitivity of the system. The change inwater level of the dryer-separator pool and the removal of the well shield plug were clearly identified.
The spatial resolution was found to be 0.5 m. GEANT simulations tuned to reproduce the observed
profiles predicted that a 2 m (cubic) fuel material should be identifiable within a week. By operating multiple systems simultaneously, or moving the system closer to the reactor, detection of 1.3 m(cubic) fuel material should also be achievable within a month.
19/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

Number of events

PTEP 2013, 073C01

PTEP 2013, 073C01

H. Fujii et al.

The constructed system operated reliably over a period of seven months, during which an interruption (due to a malfunction of the clock generator) only occurred once. The compactness of thedetector—fitting in a shipping container—is also to be noted. The demonstrated ability to investigate the inner structure of the reactor, together with its compactness and reliable operation, are very
promising features in the adoption of a muon radiography system to image the inner structure of
nuclear plants.

Acknowledgements

References[1] L. W. Alvarez et al., Science 167, 832 (1970).
[2] H. Tanaka et al., Nucl. Instrum. Methods Phys. Res., Sect. A 507, 657 (2003).
[3] K. Nagamine et al., Nucl. Instrum. Methods Phys. Res., Sect. A 356, 585 (1995).
[4] K. Nagamine et al., Proc. Jpn. Acad. B 81, 257 (2005).
[5] H. Tanaka et al., Hyperfine Interact. 138, 521 (2001).
[6] K. Morishima et al., Proc. 26th Workshop on Radiation Detectors and their Uses (KEK Proceedings
2012-8), p. 27 (2012).
[7] Plastic scintillator EJ-200 Elgin Technology.
[8] Kuraray Co., Ltd., homepage: http://www.kuraray.co.jp/en/.
[9] Hamamatsu Photonics K.K., homepage: http://www.hamamatsu.com/.
[10] A. Vacheret et al., Nucl. Instrum. Methods Phys. Res., Sect. A 623, 201 (2010).
[11] Bee-Beans Technology, MPPC Readout Unit, Catalog No. BBT-X036.nit.
[12] Bee-Beans Technology, NIM Clock distributor for BBT-X036.
[13] NEOMAX Engineering Co., Ltd., homepage: http://www.nxe.co.jp/.
[14] “Overview of Light Water Nuclear Power Stations”, edited by Nuclear Safety Research Association,(2010); http://www.rist.or.jp/atomica/data/pict/02/02030101/11.gif[15] Geant4 homepage: http://www.geant4.org/geant4/.
[16] H. Jokisch et al., Phys. Rev. D 19, 1368 (1979).
[17] S. Tsuji, Ph.D. Thesis, Okayama University (1998) [in Japanese].
[18] The fuel structure is taken retrieved from http://ousar.lib.okayama-u.ac.jp/metadata/8238

20/20

Downloaded from http://ptep.oxfordjournals.org/ at UNIV OF TSUKUBA on September 12, 2014

We would like to thank the Japan Atomic Power Company for permitting us to locate the detector systemon their site. Without their understanding, we would not have been able to accomplish the present study. We
wish to particularly acknowledge H. Okuda for informative discussions on conducting the measurement at the
JAPC whilst adhering to the Nuclear Source Material regulations. H. Yokomizo of the Japan Atomic Energy
Agency and H. Sugawara of the Okinawa Institute of Science and Technology Graduate University are also
to be acknowledged for their support and advice for conducting the study at the JAPC. The detector system
was designed and constructed in a very short period and this was realized through the excellent work by thecompanies involved and the help provided by many people at KEK; M. Tanaka for the MPPC DAQ design;
Y. Yasu for DAQ online programming; and J. Haba for budgetary management.