Monte Carlo simulation for background study of geophysical inspection with cosmic-ray muons

. Niita, K., Sato, T., Iwase, H., Nose, H., Nakashima, H. & Sihver, L., 2006. PHITS: a particle and heavy ion transport code system, Radiat. Meas., 41, 1080–1090. Nishiyama, R., Miyamoto, S. & Naganawa, N., 2014. Experimental study of source of background noise in muon radiography using emulsion film detectors, Geosci. Instrum. Methods Data Syst., 3(1), 29–39. Nishiyama, R., Miyamoto, S. & Naganawa, N., 2015. Application of Emulsion Cloud Chamber to cosmic-ray muon radiography, Radiat. Meas., 83,


Density imaging with cosmic-ray muons
Several geophysical exploration methods have been developed and employed to image the internal structure of volcanoes. An important objective for volcanologists is to put constraints on models of magma plumbing systems. From the view point of volcano monitoring, the main ambition is to detect movements of magma in magma chambers or magma pathways (conduits), for which the spatial resolution of the imaging is important because of the highly heterogeneous internal structure of volcanoes.
Recently, a novel technique named muon radiography (muography) has been developed for probing the internal density profiles of volcanoes (e.g. Tanaka et al. 2007;Ambrosi et al. 2011;Lesparre et al. 2012;Cârloganu et al. 2013). Muography is based on measurements of the absorption of the atmospheric muon flux inside a target material. The method is made possible by the fact that the energy spectrum of cosmic-ray muons and the energy dependence of muon range have been intensively studied. Thus the attenuation of the muon flux can be used to derive the column density of target mountains along muon trajectories.
The greatest advantage of muography is its high spatial resolution compared with other geophysical exploration methods. The resolution of muography can be expressed ideally as where L is the distance between the muography detector and the target, and θ is the angular resolution of the detector. Considering typical configurations (L ∼ 1 km and θ ∼ 30 mrad), the spatial resolution becomes x ∼ 30 m, better than many other geophysical methods. It is expected that muography will allow visualization of highly heterogeneous structures near the surface of volcanoes. Several types of detectors have been utilized in muography, such as scintillation detectors (e.g. Tanaka et al. 2009;Ambrosi et al. 2011;Lesparre et al. 2010Lesparre et al. , 2012Tanaka et al. 2014), gaseous detectors (e.g. Alvarez et al. 1970; Barnaföldi et al. 2012;Cârloganu et al. 2013) and nuclear emulsions (e.g. Tanaka et al. 2007;Dedenko et al. 2014;Nishiyama et al. 2014). The advantage of scintillation detectors and gas chamber detectors is that they can provide information on the arrival time of incident muons. This enables us to monitor the time variation of the muon flux and hence the sequential change of the density profile inside volcanoes. A disadvantage is that the experiment site is limited to places with well-established infrastructure, since scintillation detectors require electricity and gas chamber detectors require continuous gas circulation for operation. On the other hand, nuclear emulsion can be installed almost everywhere except for areas of high temperature (>30 • C) and it does not require electricity for operation. Another advantage of emulsion is its high resolution for determining the incident angle of muons ( θ = 5-10 mrad), which promises a better spatial resolution (eq. 1).

1040
R. Nishiyama et al. The use of muography has spread not only to volcanology but also to various fields of geoscience and civil engineering. For instance, muography has been applied to non-invasive inspection of historic architectures (Minato 1986), monitoring of furnaces (Nagamine et al. 2005), detection of seismic faults (Tanaka et al. 2011) and detection of caves (Barnaföldi et al. 2012).

Background particles
Muography observations require a precise measurement of the flux of cosmic-ray muons after passing through the target mountain (we call them 'penetrating muons' in this paper). The flux of the penetrating muons is generally very small. For instance, the muon flux for nearly horizontal muons is more than two orders of magnitude smaller than the flux of vertical muons and it decreases to two orders of magnitude smaller after passing through 1 km-thick standard rock (Groom et al. 2001). This weak flux of penetrating muons presents a difficulty for muography, since a small detector and a short exposure do not provide sufficient statistics for density estimations. To obtain sufficient statistics, one has to use a larger detector or increase the exposure time, which are often not easily achieved practically. This trade-off has been discussed in other literature (e.g. Lesparre et al. 2010).
Another difficulty in muography arises from background particles, namely particles arriving at the detector that are not penetrating muons. If the background noise exceeds the weak flux of penetrating muons from the target mountain, the muon flux may be overestimated, leading to a significant underestimation of the column density. Several research groups have reported that the flux of background particles is not negligible compared with the flux of penetrating muons (Carbone et al. 2014;Jourde et al. 2013;Nishiyama et al. 2014;Ambrosino et al. 2015). For instance, Carbone et al. (2014) has observed a remarkable difference between the observed particle flux and the expected muon flux at Mt. Etna with a scintillation detector. Specifically, the flux of particles coming from the mountain exceeds the expectation by a factor of 2-10. Ambrosino et al. (2015) also reports a similar result. Although the instrumental noise level of their detector is far below the expected muon rate, they still observed non-negligible amount of signal caused by some charged particles. This excess is attributed to contamination by background particles.
Although the origin of the background particles is not fully understood, some prior works have considered possible candidates. For instance, Jourde et al. (2013) proposed that the background noise arises from particles entering the back of the detector with upwardgoing trajectories. These particles mimic the signals of penetrating muons, since upward-going particles have trajectories identical to those of downward-going muons emerging from target mountains. They have detected these upward-going particles in several locations by dedicated time-of-flight (TOF) analysis with scintillation detectors. Although the upward-going particles are part of the source of the background particles, they do not explain all the excess in the particle flux. Nishiyama et al. (2014) proposed that the background particles are due to low-energy charged particles (E < 1 GeV) which are scattered onto detectors from random directions. They installed two emulsion detectors with different energy thresholds (0.1 and 1.0 GeV) at the foot of a small lava dome (Mt. Showa-Shinzan) and demonstrated that only the detector with the higher energy threshold yields appropriate particle flux and density values for the lava. However, the origins of these low-energy particles have not been verified, nor has their connection with the upward-going particles.
The present study aims for a better understanding of the origin of the background noise in muography via Monte Carlo simulations of particle transportation. For quantitative discussion of the background, it is required to calculate the energy spectra of major cosmic particles (not only muons) with a certain amount of systematic uncertainty. We have established a simulation framework to calculate the energy spectra of muons, electrons, gamma-rays, protons and neutrons arriving at a detector with the 3-D topography of a mountain taken into account (Section 2). The results of the simulation are shown in Section 3. The systematic uncertainty of the simulation is discussed in Section 3.2. The simulation results are compared with actual observations at Mt. Showa-Shinzan with nuclear emulsions in Section 4. In this paper, instrumental noise, such as accidental coincidence of dark noise with photo-sensors, is assumed to be sufficiently reduced and is not discussed.

S I M U L AT I O N
In this section, we describe the framework of the Monte Carlo simulation. The aim of the framework is to calculate the composition and the energy spectrum of cosmic particles arriving at a detector with the topography of a mountain taken into consideration. The calculation is performed in two steps. First, we calculate the energy spectra of major cosmic particles using the air-shower simulation code COSMOS. Second, we simulate propagation near a detector and a mountain with the GEANT4 toolkit (Agostinelli et al. 2003). Fig. 1 shows a schematic illustration of our simulation framework.

Air-shower simulation with COSMOS
COSMOS follows the development and evolution of air shower particles and describes the angular, spatial, temporal and 3-D momentum distribution of secondary particles (Roh et al. 2013). Counting the number of secondary particles enables us to collect the energy spectrum of any cosmic particle for any altitude in the atmosphere. In this study, the version of 7.641 was used (http://cosmos.n.kanagawa-u.ac.jp/).
In this work, protons and helium nuclei are injected at an altitude of 400 km. The energy distribution of the incident protons and helium nuclei are taken from the BESS flight in 1998 (Sanuki et al. 2000). The physical properties (density and composition) of the atmosphere are taken from the US Standard Atmosphere (U. S. Goverment Printing Office 1976). The spherical shape of the Earth and the atmosphere is implemented properly to calculate the flux of nearly horizontal particles. The magnetic field of the Earth is neglected in this study, although it must be implemented to discuss the longitude and latitude dependence. The absence of the geomagnetic field does not affect the general conclusion of this work, as we discuss in Section 3.2. The interaction models employed in this study are PHITS 1 (Particles and Heavy Ion Transport code System, Niita et al. 2006) and JAM (Jet AA Microscopic Transport, Nara et al. 1999) 2 for the low-energy region (E < 4 GeV), DPMJET-III (Roesler et al. 2001) for the middle-energy region (4 ≤ E < 1000 GeV) and QGSJET-II-03 (Ostapchenko 2006) for 1 PHITS is not an interaction model but a general purpose MC simulation code. In COSMOS some interaction models used in PHITS can be imported.  the high-energy region (1000 GeV ≤ E). The development (secondary production, decay, etc.) of the air-showers is traced until the kinetic energy of the secondaries drops below 50 MeV to save computation time and get better statistics. While tracing, particles passing through virtual spheres at several altitudes are recorded. The particle type, position, direction and kinetic energy are stored in the output.

COSMOS results
From the COSMOS output, we derived the energy spectrum of particle-i (i: muon, electron, gamma-ray, proton and neutron) by averaging the number of hits over the entire surface of the virtual sphere. The procedure is as follows. The number of hits is binned as a function of kinetic energy (E, index j) and zenith angle (θ , index k): where n i (···) denotes the number of hits of the ith particle, T is the equivalent exposure time and A is the geometrical acceptance for particles entering the sphere with zenith angle within θ minθ max . The radius of the Earth (R = 6.4 × 10 6 m) gives an exact value of the acceptance: A = 2π 2 R 2 (cos 2θ min k − cos 2θ max k ). α is a scaling constant and is set to 0.65 so as to fit the energy spectrum data for vertical muons (Haino et al. 2004). The difference of the normalization constant to 1 hints to the level of uncertainty affecting the simulation (∼40 per cent). Fig. 2 represents the resultant energy spectra of muons (μ ± ), electrons (e ± ), protons (p), gamma-rays (γ ) and neutrons (n) along with energy spectra reported in other literature. Although the COS-MOS results agree with the literature values in general, there are discrepancies in low-energy regions for p and n and in high-energy regions for e ± . This difference should be regarded as systematic uncertainty in the COSMOS calculation and its effect on background estimation will be discussed in Section 3.2.
For GEANT4 simulation, we produced an energy spectrum model for each particle by interpolating or extrapolating the results for 300 m above sea level (asl) (Fig. 3). The energy range of the model is {E: 1 ≤ E < 10 000 GeV} for muons and {E: 0.05 ≤ E < 500 GeV} for the other particles. The zenith dependence of the spectrum is considered by binning at intervals of cos θ = 0.05 for muons and cos θ = 0.10 for the other particles, ranging from cos θ = 0 (horizontal) to cos θ = 1 (vertical). A simple power law spectrum is used for extrapolation in the high-energy region where there are not enough statistics for fitting.

Local simulation with GEANT4
Since COSMOS cannot deal with the topography of a mountain, we use the GEANT4 toolkit to simulate particle propagation near the mountain and the detector. We constructed a virtual mountain and a virtual detector in a computational region of GEANT4 and injected particles from a substantially large hemisphere enclosing the mountain and detector, following the energy spectrum model derived from COSMOS.
The virtual mountain has a rotationally symmetric shape and is realized by a number of small prisms with horizontal dimensions x = y = 10 m (Fig. 4a). The elevation at each point (h) is given as a function of the distance to the axis (r): This shape is adjusted so that the rock thickness becomes comparable to the case of our emulsion observations at Mt. Showa-Shinzan. The virtual mountain consists of standard rock with a density of 2.00 g cm −3 . The computational region outside the terrain is filled with air (1 atm). The virtual detector has a belt-like surface and is placed surrounding the mountain at height of 65 m. The radius and height of the belt are 500 and 10 m, respectively. Thus the total area of the virtual detector is S MC = 3.1 × 10 4 m 2 . This virtual detector records the information of particles passing through it. The hemisphere, from which particles are injected, has an oblate spheroidal shape with a long axial radius of R x = R y = 600 m and a short axial radius of R z = 300 m. The size of the spheroid is adjusted so that it encloses the mountain and the detector. The particle type, position, direction and kinetic energy of incident particles are sampled based on the energy spectrum model produced by the COS-MOS simulation. The Fritiof string model (E > 10 GeV) and Bertini cascade model (E < 10 GeV) are employed as hadronic interaction models (FTFP_BERT in GEANT4 reference physics list). For electromagnetic process and multiple Coulomb scattering, the standard electromagnetic interaction model of GEANT4 was employed. To  Allkofer et al. (1985) at sea level; e + + e − : Golden et al. (1995) at 945 g cm −2 (600 m asl); p: Brooke & Wolfendale (1964) at sea level; γ : Beuermann & Wibberenz (1968) at 760 g cm −2 (2500 m asl); n: Gordon et al. (2004) at sea level (167 m asl). The solid circles denote the results of a COSMOS simulation taken at the same altitude at the experimental data for comparison.
save computation time, neutrons with kinetic energy below 30 MeV are discarded during tracing.
The rotationally symmetric mountain and detector allowed to enlarge the detector size without losing generality. The large detector acceptance significantly increased statistics of background particles with limited computation time. Specifically, we injected only 3.3 × 10 9 particles, which corresponded to the number of particles incident on the hemisphere in 10.8 s(=T MC ). However, we obtained enough statistics from the simulation, since the area of the detector (S MC ) was large. The effective exposure of the simulation The computation time was 1.6 × 10 4 hr for a single thread. With the aid of multithreading technology, the calculation was finished within one day in our computational resource. The rock thickness in R3 is 579-917 m (Fig. 7b). Fig. 4(c) shows the simulated number distribution as a function of kinetic energy of the particles when they reach the detector. In this figure, the contributions from penetrating muons and background particles are individually drawn. In both regions, the distribution of the penetrating muons shows a maximum at around 100 GeV. The penetrating muons make almost no contribution below 1 GeV, whereas the background particles (protons, electron and muons backgrounds) dominate the population in this range. We can calculate the flux of penetrating muons and background particles by integrating the number distribution over energy. The results of this calculation show that the flux of these background particles above 50 MeV in R2 is 7.8 times that of the penetrating muons, and is as much as 16.5 times in R3. This result indicates that the signals of penetrating muons would be overwhelmed by the massive flux of background particles in the case where the energy threshold of the detector is less than 1 GeV. To reduce the contamination by background particles, the optimal energy threshold should be above 1 GeV. This conclusion will be confirmed by our emulsion experiments in Section 4.

Uncertainty of simulation
In this subsection, the systematic uncertainty in the calculation of background flux is discussed.
First, we have to address the systematic uncertainty in the hadronic interaction models. Although it is very difficult to declare how the model uncertainty propagates to that of our background estimation, it can be estimated pessimistically by focusing on the discrepancy between the energy spectra derived from COSMOS and the literature values (Fig. 2). Regarding this difference as systematic uncertainty of the COSMOS simulation, the uncertainty of the background flux is estimated to be ∼40 per cent.
Second, it has to be taken into account that the magnetic field of the Earth is neglected in this calculation. There are two effects which must be considered: the overestimation of the background flux because of neglecting the geomagnetic cut-off and the underestimation of the scattered flux due to the propagation of the low-energy particles in the magnetic field. In conclusion, these effects are not so severe than the hadronic uncertainty stated above. The reasons are as follows.
(i) The geomagnetic field prevents low-energy primaries from penetrating through the magnetosphere to the atmosphere of the Earth (geomagnetic rigidity cut-off). The absence of the geomagnetic field, therefore, overestimates the flux of protons coming into the atmosphere and hence overestimates the background flux. However, the overestimation is estimated to be no more than 17 per cent, considering the vertical rigidity cut-off of Showa-Shinzan region (8 GV).
(ii) The absence of the geomagnetic field does not influence the COSMOS simulation in the top of the atmosphere because of the short length (15-30 km). The rigidity for a gyroradius of 30 km is merely ∼1 GV, assuming that the strength of the geomagnetic field is 4 × 10 −5 Tesla.
(iii) The absence of the geomagnetic field does not affect the GEANT4 simulation near the surface because of the small injection hemisphere (∼500 m). The rigidity for a gyroradius of 500 m would be the order of 10 MV. Even low-energy background particles will not be bent in the hemisphere.

Origin of background
From our simulation, the background particles can be classified in to three categories according to their origins: (i) protons, (ii) electrons and muons produced by hadronic interaction of protons and neutrons in the atmosphere and the topographic material and (iii) electrons and muons scattered in the atmosphere. In this paper, we refer to (i) and (ii) as hadronic backgrounds and we refer to (iii) as scattered backgrounds. Fig. 5 shows each component of the background particles as a function of the energy threshold for the R2 and R3 regions. In both regions, the dominant contribution is from hadronic backgrounds. The proportion of hadronic background in the total background above 50 MeV is 89 per cent and 84 per cent for R2 and R3, respectively.

Upward-going particles
Although we inject only downward-going particles (cos θ > 0) in our simulation, we find upward-going particles arriving at the detector from the rear side. The flux of the upward-going particles above 50 MeV is 5.1 × 10 −2 m −2 s −1 sr −1 (23 per cent) and 8.4 × 10 −2 m −2 s −1 sr −1 (44 per cent) for R2 and R3, respectively (Fig. 5). A zenith angle distribution of the simulated flux for the background particles is represented in Fig. 6.

Muography at Mt. Showa-Shinzan
In this section, we compare the results of our simulation with the particle flux observed with nuclear emulsions placed at Mt. Showa-Shinzan. Mt. Showa-Shinzan is a lava dome which was extruded as a parasitic cone of the Usu volcano in 1944-1945. The relative height of the peak is ∼200 m and the diameter is ∼1 km (Fig. 7a). For a comparative study, two types of emulsion detectors were installed at the foot of the mountain, 500 m west to the summit. One is a stack of four emulsion films (Quartet detector). The other is a stack of 20 emulsion films and nine 1-mm-thick lead plates, the so-called Emulsion Cloud Chamber (ECC detector;De Serio et al., 2003). The emulsion film used in this study is OPERA-type (Nakamura et al. 2006). The effective area of the two detectors is S OBS = 104 cm 2 . The exposure duration is T OBS = 1.45 × 10 7 s (168 d). The thickness of the rock existing along the muon trajectories is given in Fig. 7(b).
After chemical development of the films, the tracks in the films were scanned using the European Scanning System (Arrabito et al. 2006). The track reconstruction was conducted using the FEDRA system (Tioukov et al. 2006). We then applied our original track selection on the basis of (i) straightness and (ii) grain density (GD) of the tracks. The track selection is briefly explained in the next two subsections. Details are also given in our previous reports (Nishiyama et al. 2014(Nishiyama et al. , 2015.

Straightness cut
Stacking several emulsion films enables us to impose an energy threshold on the incident particles by analysing the straightness of the tracks (De Serio et al. 2003). The straightness is evaluated by using the deflection angle of the tracks in the adjacent films. Owing to the lead plates inserted between the films in ECC detector, this 'straightness' leads to the ECC detector gaining a higher energy threshold (∼1 GeV) than the Quartet detector (∼50 MeV). The energy dependence of the detection efficiency, calculated by a GEANT4 simulation, is presented in Fig. 8.

Grain density cut
The GD refers to the number density of silver grains produced along the tracks. It positively correlates with the ionizing power of the incident particles and hence can be used for particle identification (Toshito et al. 2004). For instance, the energy loss by ionization for particles with a kinetic energy of 100 MeV is 1.7 MeV g −1 cm 2 for electrons, 1.4 MeV g −1 cm 2 for muons and 15 MeV g −1 cm 2 for protons. Thus we can roughly distinguish thick tracks of protons from thin tracks of the others. Fig. 9(a) shows a GD distribution for selected tracks in the R2 and R3 regions of the Quartet detector. The distribution clearly has two groups: low and high GD groups. Considering the ionization power of major charged particles at near the energy threshold (∼100 MeV), the low GD group is composed of relativistic particles (protons above 1 GeV, muons and electrons) and the high GD group is composed of non-relativistic particles (protons below 1 GeV). Conversely, the Quartet detector is capable of discriminating lowenergy protons from others. In contrast to the Quartet detector, the GD distribution of the ECC detector yields a single peak because low-energy protons (<1 GeV) were already rejected based on the straightness cut. Fig. 10(a) shows the angular distribution of the tracks surviving the cuts for the low GD group of the Quartet detector, the high GD group of the Quartet detector, and the ECC detector. The particle Figure 6. Zenith angle dependence of the calculated flux of the background particles for kinetic energy above 50, 100, 200 and 500 MeV (protons, electrons and muons are added). cos θ < 0(>0) indicates downward (upward) going particle. flux was derived from the number of tracks after efficiency correction and normalization with respect to detection area and exposure time. Fig. 10(b) shows the detection efficiency for each group (see Nishiyama et al. 2014 for the Quartet detector and Nishiyama et al. Figure 8. Probability that particles pass through the emulsion detectors and are recognized as a straight track, calculated using GEANT4. This survival probability is represented for each pair of particle type and detector type. 2015 for the ECC detector). The detection efficiency becomes maximum for tracks perpendicular to the film and decreases gradually as the slope of the track (θ) increases. For the ECC detector, the efficiency is 71 per cent at θ = 0 (centre of the view), decreasing to 40 per cent at θ = 0.6. This trend is the same for the low GD tracks in the Quartet detector. The detection efficiency of the high GD tracks (protons) is higher than 85 per cent over the entire the view. Fig. 10(c) shows the particle flux for the three groups of tracks after efficiency correction. The values of the flux in each region (R1, R2 and R3) are tabulated in Table 1. For comparison, the simulation spectra convoluted with the detector response (Fig. 8) are also shown.

Observed flux
In this paragraph, we compare the observed and calculated flux for each of the three groups. The flux derived from the ECC detector agrees well with the MC expected flux (first row in Table 1). Although there is a slight deviation beyond the observational error in R3, this difference can be attributed to the fact that the mountain realized in the MC computational region is quite simplified and does not reproduce the real topography of Mt. Showa-Shinzan.    The flux derived from the high GD group of the Quartet detector should be compared to the MC expectation of the proton flux (second row in Table 1). The values of the observed and expected flux agree well within the systematic uncertainty of the simulation (∼40 per cent, see Section 3.2). It is noticeable that a substantial amount of protons come from the direction of the mountain. The particle flux derived from the low GD group of the Quartet detector should be compared with the sum of the MC expectation for muons and electrons (third row in Table 1). In the free sky (R1), both the observed and calculated fluxes are higher than the flux from the ECC detector. This is because low-energy cosmic electrons (>50 MeV) contribute approximately 10 per cent of the observed 1048 R. Nishiyama et al. flux. It is clear from the results for R2 and R3 regions that background particles have a significantly higher flux than the penetrating muons.
The particle flux derived from the Quartet low GD group exceeds significantly the flux derived from the ECC detector. The ratio Quartet low GD/ECC is 5.2 and 7.3 for R2 and R3, respectively. These values agree well with the MC expectations, 4.0 and 7.4 for R2 and R3. If we assume that all the observed tracks in the Quartet detector are due to penetrating muons, the muon flux is significantly overestimated and hence the target density is significantly underestimated. Fig. 11 shows the estimated density of Mt. Showa-Shinzan from the Quartet and ECC detectors. The density values derived from the ECC detector (1.98-2.39 g cm −3 ) are consistent with the bulk density of typical volcanic rocks. On the other hand, the density values derived from the Quartet detector are significantly lower, lower even than the water density (1 g cm −3 ) in the lower elevation regions. This fact strongly suggests that the ECC detector collects only penetrating muons and the Quartet detector is affected by contamination from low-energy (<1 GeV) charged particles other than penetrating muons.

D I S C U S S I O N A N D C O N C L U S I O N
From the numerical and experimental studies, it is concluded that the origin of the background noise is low-energy charged particles such as electrons, muons and protons, in case of muographic measurements performed with muon trackers without particle identification capabilities and with energy thresholds below 1 GeV. In this section, this result is compared with other prior works. Carbone et al. (2014) have observed an excess of particle flux of a factor of 2-10 in observations at Mt. Etna. Our simulation and observation agree with this work in this point. However, it should be kept in mind that the Etna observation was performed at ∼3000 m asl and our simulation assumes nearly sea level observation. The flux of electrons at 3000 m asl would be one order of magnitude higher than that at sea level. In addition, the flux of protons and neutrons depends more strongly on the altitude. Further consideration is required to discuss the altitude dependence of background flux. In addition, the background flux depends not only on altitude but also on the longitude and latitude of the detector site through a rigidity cut on primary cosmic rays, but this effect was not considered in our simulation. The geomagnetic effect must be properly implemented in the simulation in future studies for a more quantitative discussion. Jourde et al. (2013) experimentally demonstrated the existence of upward-going particles which enter the detector from the rear side. Our simulation also verified the existence of upward-going particles. However, it must be mentioned that the backgrounds consist not only of upward-going particles but also of downward-going particles with significant abundance. As seen in Section 3.4, the upward-going particles amount for only 44 per cent of the total flux of background particles in the cos θ range of [0: 0.15], and 23 per cent in [0.15: 0.25]. The flux of upward-going particles decreases with an increase of elevation angle. This result is quite consistent with Jourde et al. (2013) as that work reports that the upward-going particles are observed mainly at low elevation angles.
This work measures the bulk density of Showa-Shinzan lava to range between 1.98-2.39 g cm −3 . The associated error (1σ ) is 0.07-0.23 g cm −3 . By contrast, Tanaka et al. (2007) performed muography observation at the same target with emulsion detectors and obtained 2.71-2.91 g cm −3 with a 1σ error of 0.17 g cm −3 at each data point. The Showa-Shinzan lava is classified as dacite from the SiO 2 content of rock samples. The bulk density of the rock samples is 2.32 g cm −3 (Nemoto et al. 1957). The prior work exceeds the sampled density significantly. On the other hand, the density values estimated from the ECC detector of the present work agree well with the sampled density. This shows that the accuracy of muography has been improved with the aid of background reduction with the ECC detector. The detector used in the prior work was a stack of four emulsion films. Thus this observation would have been affected by background noise. To avoid underestimation of the density as seen in the case of the Quartet detector, this prior work performed background subtraction (Okubo & Tanaka 2012). To be precise, the flux of background particles was estimated from the flux of nearly horizontal directions and was subtracted from the observed flux in the other directions. This subtraction process is based on an assumption that the flux of background particles is isotropic. However, this assumption is not correct. The background flux calculated by simulation and measured by emulsions depends strongly on incident angle.
In conclusion, the origin of the background noise of muography has been confirmed to be (i) protons, (ii) muons and electrons produced by hadronic interaction of protons and neutrons in the atmosphere and the topographic material, and (iii) muons and electrons which are scattered in the atmosphere. The flux of these particles exceeds the weak flux of penetrating muons at low energy (<1 GeV). Thus the only way to achieve accurate muography is to impose an energy threshold on the particles entering detectors. It is not sufficient to reject upward-going particles by TOF analysis, since the background contains downward-going components as well as upward-going particles. It is still not satisfactory to perform particle identification, since the backgrounds contain low-energy muons produced by hadronic interaction of protons and neutrons. The present work confirms that a desirable energy threshold is 1 GeV when the observation is performed at nearly sea level and the target thickness is approximately 1 km. The result will be of importance for detector design in future works, not only for emulsion detectors but also for other types of detector. At the same time, previous muography works should be reviewed from the viewpoint of whether or not they have been affected by background particles.

A C K N O W L E D G E M E N T S
For this study, we used the computer systems of the Earthquake and Volcano Information Center of the Earthquake Research Institute, University of Tokyo. This work greatly benefited from discussions and encouragement from H. K. M. Tanaka and S. Okubo (The University of Tokyo). RN has been supported by a fellowship and grant from JSPS (DC2-25-9317) during the essential part of this study. We would like to thank N. Mitsuhiro (F-lab, Nagoya University), C. Bozza (Salerno University) and V. Tioukov (INFN, Napoli) for their support on the emulsion experiments. We would like to thank H. Oshima and T. Maekawa (Hokkaido University) for their support with the field work at Mt. Showa-Shinzan. We thank two anonymous reviewers for valuable comments to improve the manuscript.