Fast Muon Simulation in the JUNO Central Detector

2.3 MB·

by  Tao Lin

Submitted to ‘Chinese Physics C’
Fast Muon Simulation in the JUNO Central Detector *
Lin Tao 1;1 Deng Zi-Yan 1 Li Wei-Dong 1 Cao Guo-Fu 1
You Zheng-Yun 2 Li Xin-Ying 3
1 InstituteofHighEnergyPhysics,ChineseAcademyofSciences,Beijing100049,China
2 SchoolofPhysics,SunYat-senUniversity,Guangzhou510275,China
3 SchoolofPhysics,NankaiUniversity,Tianjin300071,China
6
1
0 Abstract: The Jiangmen Underground Neutrino Observatory (JUNO) is a multi-purpose neutrino experiment
2 designed to measure the neutrino mass hierarchy using a central detector (CD), which contains 20 kton liquid
scintillator (LS) surrounded by about 17,000 photomultiplier tubes (PMTs). Due to the large fiducial volume and
n
a huge number of PMTs, the simulation of a muon particle passing through the CD with the Geant4 toolkit becomes
J an extremely computation-intensive task. This paper presents a fast simulation implementation using a so-called

0 voxel method: for scintillation photons generated in a certain LS voxel, the PMT’s response is produced beforehand
3 with Geant4 and then introduced into the simulation at runtime. This parameterisation method successfully speeds

up the most CPU consuming process, the optical photon’s propagation in the LS, by a factor of 50. In the paper,
]
t the comparison of physics performance between fast and full simulation is also given.
e
d
Key words: JUNO, Central Detector, Fast Simulation, Geant4
-
s
PACS: 29.40.Mc, 29.85.Fj
n
i
.
s 1 Introduction 20kt liquid scintillator (LS) and surrounded by about
c
i 17,000 20-inch photomultiplier tubes (PMTs). An en-
ys The Jiangmen Underground Neutrino Observatory ergy resolution of 3%/(cid:112)E(MeV) is expected. Around
h (JUNO)[1, 2] is a multiple purpose neutrino experiment the CD, there is a water pool to shield radioactivities.
p todetermineneutrinomasshierarchyandpreciselymea- The PMTs in the water pool are used to detect the
[ sure oscillation parameters. It is being built in Jiang- Cerenkov lights yielded by cosmic ray muons. On top

1 men in Southern China and is about 53km away from of the water pool, there is a top tracker, made of plastic
v Yangjiang and Taishan nuclear power plants, respec- scintillators, to identify muon tracks.
6 tively. To meet the requirement of high energy resolution
5
for the CD, reliable and flexible Monte Carlo (MC) sim-
0 Top Tracker
ulation software is a necessity, especially for optimizing
0
0 Water detector parameters at the design stage. The simulation
Pool
. software has been developed based on the Geant4 [3]
2
0 20 kt LS toolkit, which consists of a number of packages to man-
6 age event generators, geometry and materials, physics
1 processes, tracking and user interfaces. Both optical pa-
: Central Detector
v rametersandphysicsprocesseshavebeentunedbasedon
i 17,000 the experimental data obtained by the Daya Bay Neu-
X
20-inch PMTs trino Experiment[4]. The simulation software is highly
r
integrated with the underlying software framework of
a
SNiPER[5] so that it becomes an important component
in the full data processing chain.
Fig. 1: (color online) Schematic view of the JUNO de- 2 Data processing
tector
The schematic view of the JUNO detector is shown in The raw data recorded by the JUNO detector need
Fig.1. Themostinnerpartiscalledthecentraldetector tobeprocessedofflineinordertobeconvertedtorecon-
(CD). The CD is basically an acrylic sphere filled with structeddata,whicharesuitableforphysicsanalysis. In
∗ThisworkissupportedbytheStrategicPriorityResearchProgramoftheChineseAcademyofSciences(GrantNo. XDA10010900),
NationalNaturalScienceFoundationofChina(GrantNo. 11405279,11575224).
1)E-mail:[email protected]
1
Submitted to ‘Chinese Physics C’
the current phase, the raw data will be produced by the tionphotonsisgivenbytheBirks’Law[9]. Thenumber
MC simulation and then are processed in the same way of scintillation photons, dS, as a function of the visible
as the real data in the future. Following the scheme of energy dE is given by
vis
SNiPER, various data processing steps shown in Fig. 2
dE dE
are implemented as corresponding type of algorithms, dS=A·dE =A ,δ= [MeVg−1cm2],
softwarecomponentsundertakingcertainamountofcal- vis 1+C δ+C δ2 ρdx
1 2
(1)
culation task. Event data pass between different algo-
where dE is the energy loss, A is the light yield and C
rithms via a data buffer. The event data cached in the 1
and C are Birks constants.
data buffer can be written to ROOT [6] files for persis- 2
After an optical photon is emitted in the LS, it will
tent storage. The detector parameters stored in GDML
travel in the LS and there is a certain possibility for it
[7] files can be accessed by algorithms through the de-
to reach a PMT after passing through the acrylic sphere
tector description service in a uniform way.
and buffering water. During the photon’s propagation,
it might be absorbed or scattered. At the liquid scintil-
lator and acrylic sphere boundary, it may be refracted,
orreflected,orevenatotalinternalreflectionmayoccur.
If it reaches the photo-cathode of a PMT, a photoelec-
tron may be generated. The timing information of the
photoelectron can be described by Eq. 2, where T is
start
theinitialtime,∆T isthegenerationtimeofthescin-
gen
tillation photon, ∆T is called hit time, which includes
Fig. 2: Work flow for data processing hit
theopticalphoton’spropagationtimeandthephotoelec-
As shown in Fig. 2, the physics generator generates tron’s generation time.
simulated physics events such as inverse β decay (IBD)
events,cosmicraymuonevents,radioactivityeventsand T =Tstart+∆Tgen+∆Thit. (2)
soon. ThegeneratedGenEventobjects,whichareinthe
format of HepMC [8], mainly include kinematics infor- The MC simulation of muon particles is important
mation of the generated final-state particles. Then de- for physics studies, but the full simulation with Geant4
tector simulation loads GenEvent objects from the data is extremely CPU consuming. To simulate a muon with
buffer and starts tracking with Geant4. The informa- typical energy passing through the central detector, at
tion of simulated hits such as charge and time informa- least 50 minutes are required using a modern CPU core.
tion are saved into SimEvent objects. After inputting So the generation of cosmic ray muons with full simula-
SimEvent objects, electronics simulation performs digi- tion becomes a big challenge. The most time consuming
tization and generates ElecEvent objects containing in- part is the propagation of optical photons due to the
formation of waveforms with predefined sampling rate. huge number of photons (at the level of 107).
Then the PmtRec algorithm processes the waveforms
and converts them into CalibEvent objects using cali- 4 Fast muon simulation
bration constants. Finally, the event reconstruction al-
gorithm reads in CalibEvent objects, performs event re- The scintillation photons make up the majority of
constructing and stores the produced RecEvent objects opticalphotonsgeneratedbythemuonspassingthrough
in the data buffer. theCD.Tospeeduppropagationprocessofscintillation
photons, a fast simulation employing a so-called voxel
3 Challenge of muon simulation method has been implemented. Instead of the full simu-
lation, it models the response of the PMTs by sampling
The muon induced background is one of the main the response distributions prepared in advance with the
backgrounds in the JUNO experiment. The mean en- Geant4 simulation.
ergy of muons that penetrate the overburden is about
4.1 Voxel method
215GeV and they reach the JUNO detector at the rate
of about 3Hz [2]. In the voxel method, the volume of the whole liq-
When a muon travels through the CD, energy is de- uid scintillator sphere is treated as a regular grid and
posited in the LS and both scintillation and Cerenkov a voxel represents the volume element in the grid. By
photons are emitted: the formers are emitted isotropi- dividing the sphere into voxels, a muon particle encoun-
cally, while the latters are produced with a fixed angle ters a serial of voxels along its path when it enters the
withrespecttothemuontrack. TheratioCerenkovover CD. In Fig. 3, the scintillation photons generated in the
scintillation photons is only 5%. The yield for scintilla- shaded voxel (fired voxel) reach the PMT and generate
2
Submitted to ‘Chinese Physics C’
detectorresponse. Thelocationrelationshipbetweenthe larger than the normal hit time.
fired voxel and the PMT can uniquely be determined by
a pair of variables (R, θ), where R is the radius of the
voxel and θ is the spatial angle between the voxel and
the PMT. If there is a certain amount visible energy of
Evis ina voxel, theresponse ofthe PMTswiththe same ns300
(R, θ) pair value follows the same distributions. It is me/250
obviousthattheresponsedistributionscanbeeasilyob- hit ti125000
tained by running the full simulation. So the core of the
100
voxelmethodistobuildtheconnectionbetweentheE
vis 50
in a voxel and the response of the PMTs. 0
3
muon Water θ/rad2.5 21.5 10.5 0 0 1000 2000 3000 4000 5R030/0m3
vertex LS
R Fig. 4: Profile of mean hit time for scintillation photons
PMT within a bin in the R-θ parameter space.
θ
center
At the stage of using the prepared response distribu-
Fig. 3: Illustration of mapping the visible energy within tionstodothefastmuonsimulation,theprocedureisas
a voxel to the response of the PMTs follows: foreverystepofmuontrackinginLS,thevisible
energy, E , is calculated according to the Birks’ Law,
Atthestageofpreparingresponsedistributions,pho- vis
which is the same as that in scintillation process. If the
tonsareproduceduniformlyinallLSvoxelswiththefull
E value is none-zero, the voxel that the tracking step
simulation. Since photons with different wavelength be- vis
encounters is identified by the step’s position and the
have differently in the detector, the emission spectrum
two associated nPE and hit time histograms are loaded.
of the LS is used to calculate the wavelength of the pho-
Then E is broken into an integer part (floor of E ,
tons. Foreachvoxel,thetotalnumberofgeneratedpho- vis vis
(cid:98)E (cid:99)) and a fractional part. Both are used to calculate
tons is equivalent to the light yield of 1MeV visible en- vis
nPE, as Eq. 3 shows, where N is the integer part and
ergy, which makes it possible for response distributions int
N isthefractionalpart. Fortheintegerpart,thenPE
to contain sufficient information for the fast simulation frac
histogramissampled(cid:98)E (cid:99)times,asshowninEq.4. For
to simulate a physics event. Photons received by each vis
the fractional part, the nPE histogram is sampled only
PMTandtheassociatedtiminginformationarethemost
once to get N(cid:98)Evis(cid:99). N follows the binomial distribu-
importantdatatoreflectopticalphoton’stransportation int frac
tion N ∼B(N(cid:98)Evis(cid:99),p), where the fractional part is
processintheLS.Sothenumberofphotoelectrons(nPE) frac int
used to calculate the probability for the PMT to be hit
andhittimecanberegardedasthemajorparametersof
or not.
PMT response.
Forsimplicity,anR-θparameterspaceisusedtorep- N=N +N (3)
int frac
resent the set of possible combinations of (R, θ) values,
in which R ranges from 0 to the inside radius of the
acrylic sphere and θ is between 0 to π. Both R and θ
are divided into a serial of intervals, respectively. Each j<(cid:88)(cid:98)Evis(cid:99)
N = Nj (4)
2-dimensional bin in parameter space is associated with int int
a number of voxels and PMTs. For a specific bin, two j=0
histograms, one is for nPE distribution and the other
If the PMT’s nPE is greater than zero, its hit time is
for hit time distribution, are generated by the full sim-
obtained by sampling the associated hit time histogram.
ulation. After traversing all the bins in the parameter
space, a complete set of response distributions are ob-
4.2 Performance measurements
tained. Fig. 4 shows a profile of mean hit time for scin-
tillationphotonswithinabinintheR-θparameterspace, The data sample used in the physics performance
in which the x-axis represents R3, the y-axis represents studies are single γs with the energy of 1MeV, which
θ andthez-axisindicatesthemeanhittime. Thisfigure aregeneratedalongverticalzaxiswithastepof0.1min
showsthatvertexattheedgeofdetectorisinfluencedby the LS. To guarantee sufficient statistics, 10,000 events
the total internal reflection, so that the mean hit time is are generated at each position.
3
Submitted to ‘Chinese Physics C’
107
1400
106
full sim 1200 full sim
105
1000
fast sim fast sim
es 104 es 800
Entri 103 Entri 600
102 400
10 200
1 0
5 10 15 200 250 300 350 400 450 500
nPE hit time/ns
(a) nPE of PMT Fig. 6: Hit time of the PMTs in the region of 175◦ to
176◦ with γs generated at z=17.3m.
107
106 Fig. 7 shows the average total number of photoelec-
full sim trons per event (totalPE) as a function of γ’s generation
105
position z. The difference between fast and full simula-
fast sim
s 104 tion is less than 0.4%, which again demonstrates a good
e
Entri 103 consistency between them. The totalPE increases with
theγ’sgenerationpositiongettingclosertototheacrylic
102
sphere. After reaching the turning point, the totalPE
10 drops down mainly because of the total internal reflec-
tion and energy leak. So the fast simulation can model
1
0 500 1000 1500 2000 2500
the change of totalPE with the vertex in a precise way.
hitTime/ns
(b) Hit time of PMT
1450
1400
Fig. 5: Comparisons of nPE and hit time between fast full sim
and full simulation using γs generated at z=16m. 1350
fast sim
E
alP 1300
ot
t
1250
1200
1150
0 2 4 6 8 10 12 14 16 18
As shown in Fig. 5, the distributions of nPE and hit
z/m
timearecomparedrespectivelybetweenthefastandfull
simulation using γs generated at z=16m. It is obvious Fig. 7: Average totalPE as a function of γ’s generation
thattheagreementisgoodforbothnPEandhittimeat position z.
the PMT level. The optical photons received by a PMT
should have two components: photons directly from the When a muon particle passes through the CD, the
interaction position and reflected photons. The closer number of scintillation photons it generates is propor-
to the acrylic sphere the γs are generated, the more op- tional to its visible energy in the liquid scintillator. So
tical photons are expected to be bounced back at the the more energy the muon particle deposits, the more
barrier due to the reflection effect. For γs generated at CPU time is required to simulate the generation and
z=17.3m, Fig. 6 shows the hit time distribution of some transportation process of scintillation photons. To make
selectedPMTswhoseθvariesfrom175◦to176◦. Accord- athoroughtimingmeasurement,abroadbandofenergy
ing to MC truth information of the full simulation, the pointscovering10GeV,100GeV,215GeV,500GeVand
first peak is caused by the direct photons, however, the 1TeV is chosen for muon particle’s generation. Ahead
second one is mainly contributed by the reflected pho- of timing measurements, a full simulation job is run and
tons. So the fast simulation can reflect effectively the duringjobexecutioninformationateachMCsimulation
transportation of optical photons in the LS. steparecollectedandstoredinadatafile. Theprepared
4
Submitted to ‘Chinese Physics C’
datafilewillbeusedbythesubsequentfastsimulationso uid scintillator. Although the voxel method has been
thatitcanbecomparedwiththefullsimulationatevent developed to speed up the muon simulation in the cen-
level. The measurement of system performance is done tral detector, it can also be extended to simulate other
on a blade server with the CPU Intel(cid:13)R Xeon(cid:13)R E5-2680 types of events such as IBD because the parameteriza-
v3 @ 2.50GHz. In order to eliminate interferences, the tion approach is related to neither event types nor event
CPU is exclusively used by the timing measurement job energy.
and other irrelevant applications such as system moni- The symmetry of detector geometry greatly reduces
toring are all suspended. the number of sampling histograms used by the voxel
method. But there are many factors that might break
9 this kind of symmetry. For example supporting sticks
around the acrylic sphere will block detection of optical
8
min 7 photons. The voxel method can handle this problem by
me/ 6 E5 2680 v3 increasing the number of sampling histograms.
Ti @2.50GHz Itiseasytouseparallelcomputingtechniquesuchas
n 5
o CUDA [10] to further accelerate execution of the voxel
ati 4
ul method. Unlike the full simulation application running
m 3
Si on GPU such as Chroma [11], the parallel implementa-
ast 2 tion of voxel method can only execute histogram sam-
F
1 pling with GPUs, which avoids the complexity of geom-
50 100 150 200 250 300 etry translation.
Full Simulation Time/min
Fig. 8: The fast simulation time versus the full simula- 5 Conclusions
tiontimeformuonparticleswiththeirmomentavarying
between 10GeV and 1TeV.
Thefastmuonsimulationusingthevoxelmethodhas
Fig. 8 is a scatter plot of execution time elapsed on been implemented in the JUNO’s offline software frame-
generation and transportation of scintillation photons work. The applied method is to generate the response
generated by muons passing through the central detec- ofthePMTsbeforehand,forphotonsproducedindiffer-
tor. Eachcrossrepresentsthefastsimulationtimeversus ent voxels in the CD with Geant4 and then to introduce
the full simulation time for a muon particle. The result the response into the fast simulation by sampling the
shows that the elapsed simulation time scales linearly response histograms at runtime.
with the muon particle’s visible energy and the fast sim- The timing measurement shows that the fast simula-
ulationachievesaspeedupratioofabout50overthefull tionhasobtainedaspeedupratioof50inCPUexecution
simulation. time compared to the full simulation with Geant4. The
physics performance has also been examined and valida-
4.3 Discussions
tion results show that the agreement between the fast
The JUNO CD simulation with Geant4 is largely and full simulation is good and no significant discrepan-
dominated by optical photons’ transportation in the liq- cies have been found.
References 81-86(1997)
7 R.Chytracek,J.McCormick,W.PokorskiandG.Santin,IEEE
1 Z.Djurcicetal(JUNOCollaboration),arXiv:1508.07166 Trans.Nucl.Sci.,53: 2892(2006)
2 F.P. An et al (JUNO Collaboration), arXiv:1507.05613, ac- 8 M. Dobbs and J. B. Hansen, Comput. Phys. Commun., 134:
ceptedforpublicationinJ.Phys.G. 41-46(2001)
3 S. Agostinelli et al (GEANT4 Collaboration), Nucl. Instrum. 9 J.B.Birks,Proc.Phys.Soc.A,64: 874(1951)
Meth.A,506: 250-303(2003) 10 https://docs.nvidia.com/cuda/cuda-c-programming-guide/,
4 F.P.Anetal(DayaBayCollaboration),Phys.Rev.Lett.,108: retrieved10thMarch2015
171803(2012) 11 http://chroma.bitbucket.org/ downloads/chroma.pdf, re-
5 J.H.Zouetal,J.Phys.Conf.Ser.,664(7): 072053(2015) trieved12thDecember2015
6 R. Brun and F. Rademakers, Nucl. Instrum. Meth. A, 389:
5