- Email: [email protected]
Study of Track Reconstruction for DAMPE
CHINESE ASTRONOMY AND ASTROPHYSICS Chinese Astronomy and Astrophysics 41 (2017) 455–470
Study of Track Reconstruction for DAMPE† LU Tong-suo1,2,3
LEI Shi-jun1,2
CHANG Jin1,2 1 2
ZANG Jing-jing1,2
WU Jian1,2
Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008
Key Laboratory of Dark Matter and Space Astronomy, Chinese Academy of Sciences, Nanjing 210008 3
University of Chinese Academy of Sciences, Beijing 100049
Abstract The Dark Matter Particle Explorer (DAMPE) is aimed to study the existence and distribution of dark matter via the observation of high-energy particles in space with a large energy bandwidth, high energy resolution, and high spatial resolution. The track reconstruction is to restore the positions and angles of the incident particles using the multiple observations of different channels at different positions, and its accuracy determines the angular resolution of the explorer. The track reconstruction is mainly based on the observations of two sub-detectors, namely, the Silicon Tracker (STK) and the BGO (Bi4 Ge3 O12 ) calorimeter. In accordance with the design and structure of the two sub-detectors, and using the data collected during the beam tests and ground tests of cosmic rays, we discuss in detail the method of track reconstruction for the DAMPE, which includes mainly three basic procedures: the selection of track hits, the fitting of track hits, and the judgement of the optimal track. Since an energetic particle most probably leaves multiple hits in different layers of the STK and BGO crystals, we first provide a method to obtain a rough track in the BGO calorimeter by the centroid method, and hereby to constrain the track hits in the STK. Then for the selected one group of possible track hits in the STK, we apply two different algorithms, the Kalman filter and the least square linear fitting, to fit these track hits. The consistency of the results obtained independently via the two algorithms confirms the validity of our track reconstruction results. Finally, several criteria for picking out the most possible track among all the tracks found in the reconstruction by combining the results †
Supported by National Natural Science Foundation (11303105, 11303106, 11303107)
A translation of Acta Astron. Sin. Vol. 57, No. 3, pp. 353–365, 2016
Received 2015–10–16; revised version 2015–11–16
[email protected]
[email protected]
0275-1062/17/$-see front matter © 2017 Elsevier B.V. All rights reserved. doi:10.1016/j.chinastron.2017.08.012
456
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
of the BGO calorimeter and STK are discussed. Using the track reconstruction method proposed in this article and the beam test data, we confirm that the angular resolution of the DAMPE satisfies its design requirement. Key words instrumentation: detectors—planets and satellites: detection— gamma rays: ISM— cosmology: dark matter—methods: data analysis 1.
INTRODUCTION
Dark matter is a puzzle in astronomy and the development of modern astrophysics, now the detection of dark matter is during the period of rapid development. In China, aiming at the dark matter, the Dark Matter Particle Explorer (DAMPE) was launched into space on 17th December 2015, it adopts a 3-dimensional total-absorption calorimeter of large dynamic range, and has a series of advantages, such as the large energy range, high energy resolution, and the strong ability of back ground suppression, etc. According to the requirement of its physical design and the result of simulations, the DAMPE is composed of 4 sub-detectors, from top to bottom, they are the Plastic Scintillator Detector (PSD), Silicon Tracker (STK), BGO calorimeter, and Neutron Detector (NUD), as shown by Fig.1. The effective detection area of the PSD subsystem is 820 mm×820 mm, the whole subsystem consists of two layers of totally 82 plastic scintillator modules perpendicular to each other along the X and Y directions, in which the size of 78 plastic scintillator modules is 884 mm×28 mm×10 mm, and the size of other four modules is 884 mm×25×10 mm. On the both ends of the PSD unit, photomultipliers are adopted to convert the optical signals into electric signals for further treatment. The STK has 6 major layers, each major layer consists respectively of two layers of silicon microstrip detectors perpendicular to each other along the X and Y directions, mainly for measuring the directions and tracks of incident particles. A piece of tungsten plate of 1 mm thickness is placed respectively upon the second, third, and fourth major layers, to convert Gamma photons into positive and negative electrons. Considered that after the photon conversion, it is necessary to measure the tracks of the produced secondary particles, in order to distinguish the charged particles, especially distinguish the Gamma photons and electrons[1] , under the tungsten plate there must be multiple layers of silicon microstrip detectors[2] . AS shown by Fig.2, the main body of the BGO calorimeter is made of 308 pieces of BGO crystals, its detection area is 60 cm×60 cm. The BGO calorimeter is composed of totally 7 major layers, each major layer consists of the X and Y two minor layers. And each minor layer has 22 detector units composed of 22 pieces of crystals. The size of each detector unit (BGO crystal) is 2.5 cm×2.5 cm×60 cm. On the two ends of the detector unit, a R5610A-01 photomultiplier (PMT) is adopted to convert the optical signals into electric signals[3−6] . The NUD is positioned on the bottom of the DAMPE, it adopts the B-doped plastic
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
457
scintillator detector with a thickness of 1 cm (BC454, made by the Saint Gobian Company). The detector is composed of the plastic scintillator with a profile size of 693 mm×693 mm, which is divided into 4 independent squares, and one angle of each square is cut off for coupling with one PMT readout to form a complete detection plane. It is designed for measuring the secondary particles produced by the interactions of the hadrons in cosmic rays (mainly protons) with the matter on the upper layer of the NUD, according to the energy precipitation of these protons in the detector, the types of incident particles can be judged, and in cooperation with the BGO calorimeter, the protons and electrons can be further distinguished[7,8] .
Fig. 1 A sketch of 4 sub-detectors
2.
Fig. 2 The schematic diagram of the BGO calorimeter
PROCEDURES OF TRACK RECONSTRUCTION
Considered the geometric configuration and accuracy of the detector, in order to obtain a mathematic model to describe approximately the particles, we adopt the track fitting method shown in Fig.3. 3.
TRACK RECONSTRUCTION OF PARTICLES IN THE BGO CALORIMETER
The track reconstruction of particles in the BGO calorimeter is based on the following consideration: when the incident particles pass through the calorimeter, because of the interactions (mainly the electromagnetic interaction and strong interaction) of particles with the medium atoms, the particles will deposit energies in the calorimeter, according to the shower shape, we can judge approximately the track direction. The practical realization is as follows: at first, to select a BGO layer, by setting an energy threshold to get rid of the background noise, according to multiple statistical results, we preliminarily set the threshold
458
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
as 2%, only when the energy of a BGO layer exceeds 2% total energy, can this layer be taken as the reference layer to join in the track reconstruction in the BGO calorimeter. Secondly, to select the BGO crystals in each BGO layer to participate in the track reconstruction, we select the crystal with the greatest energy deposition and the two crystals in its left and right sides, totally 3 crystals, to calculate the centroid of energy deposition of this layer, the details are shown as Eq.(1). If in a layer the crystal with the greatest energy deposition is positioned on the edge, namely there is no crystal in its left or right side, in order to reduce errors, we select only this crystal with the greatest energy deposition to join in the track reconstruction; if in a layer, the crystal with the greatest energy deposition is not positioned on the edge, but the energies of the crystals in its right and left sides have a large difference, namely the energy of one side is negligible in respect to the other side, we still take the two crystals with the higher energies to participate the reconstruction, because the situation that a particle penetrates through simultaneously two neighboring crystals and gets out from the underneath is possible. Then, to make linear fitting: the energy-weighted linear fitting, the energy used for weighting is the energy summation of all the crystals participating in the reconstruction. Emax−1 xmax−1 + Emax xmax + Emax+1 xmax+1 = xc , Emax−1 + Emax + Emax+1
(1)
In Eq.(1), xc , xmax , and Emax indicate respectively the position of energy centroid of a layer in the BGO calorimeter, the position and energy of the crystal with the greatest energy in this layer. xmax−1 , xmax+1 , Emax−1 , and Emax+1 express respectively the positions and energies of the crystals in the left and right of xmax . Fig.4 shows the result of track reconstruction of a vertically incident proton of 400 GeV. This figure is the projective diagram on the XZ plane, the several lines of dots in the right side are the energy centroid positions of the crystals actually being hit in the calorimeter, the length of each transverse error bar indicates the magnitude of energy deposition, the straight line is the fitted particle track. In this figure the 6 points (X <0) belong to the cluster in the STK. 4. 4.1
TRACK RECONSTRUCTION OF PARTICLES IN THE STK
The Selection of Cluster
As the size of the BGO crystal (2.5 cm×2.5 cm×60 cm) is relatively large, so some errors exist unavoidably in the process of track reconstruction, this is caused by the system itself, we can use this error range to constrain the selection of cluster in the STK. Fig.5 shows the distribution of particle track direction reconstructed by using the BGO calorimeter information, which can be taken as a reference for the particle track reconstruction in the STK. The upper left panel displays the distribution of the angle θx between the projection of the particle track in the BGO calorimeter on the XZ plane and the Z axis,
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
459
Set the thresholdˈfind out the position of energy center
Fitting the track in BGO calorimeter
Selecting the clusters
Selecting the candidate tracks
Find out the best track
Kalman filter
The track Fig. 3 Flow chart for reconstructing the incident particle trajectory
300
X/mm
200 100 0 100 200 200
100
0
100 Z/mm
200
300
400
Fig. 4 The result of track reconstruction of a proton vertically incident to the top surface of the DAMPE with a momentum of 400 GeV
460
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
as a function of event accumulation; the lower right panel expresses the distribution of the angle θy between the projection of the particle track in the BGO calorimeter on the Y Z plane and the Z axis, as a function of event accumulation; in the lower left panel the abscissa expresses the angle θx between the projection of the particle track in the BGO calorimeter on the XZ plane and the Z axis, the ordinate expresses the angle θy between the projection of the particle track in the BGO calorimeter on the Y Z plane and the Z axis. The upper right panel shows the distribution of the angle θ between the particle track and the Z axis as a function of event accumulation. According to the calculation, in order to reduce the errors, the deviation of the slope (inclination) of the reconstructed track in the STK from that in the BGO calorimeter should be within ±5◦ , from Fig.5 we can find that most hit points are distributed in a 5◦ range, and the intercept is taken to be ±25 mm, namely the search of particle track in the STK is restricted in this range. We assume that the equation of projection of the reconstructed track in the BGO calorimeter on the XZ plane is l c : z = ax + b, as shown by Fig.6. Then, the equations of the straight lines which restrict the cluster search range are: l2 : z = tan(arctan a + 5)x + 25 + b ,
(2)
l1 : z = tan(arctan a − 5)x − 25 + b ,
5 4 3 2 1 0 1 2 3 4 5 4 3 2 1
1000 Entries
800 600 400 200 0
1 2 θx /(°)
3
4
5
6
0
Entries
1600 1400 1200 1000 800 600 400 200 0 4 3 2 1
θy /(°)
Entries
in which, a and b are respectively the slope (inclination) and intercept of the equation of projection.
0
1 2 θx /(°)
3
4
5
6
2
4
6
8
1800 1600 1400 1200 1000 800 600 400 200 0 5 4 3 2 1
Fig. 5 The track reconstruction in the BGO calorimeter
10 12 14 16 18 20 θ /(°)
0 1 θy /(°)
2
3
4
5
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
4.2
461
Track Search
Track search is a kind of pattern recognition, belonging to the scope of artificial intelligence[9,10] , its task is to divide a group of experimentally measured points into two types: the first type can be divided into many subtypes, the points in each subtype are all caused by one and the same particle, namely, corresponding to one track; the second type indicates those points which can not be attributed to any particle under sufficient confidence level, namely the noise. At present, the method of track recognition is of two kinds, namely the localized method and overall method. In the localized method, each time only one candidate track is selected, typically it starts from a few points (the initialization of candidate track), then it makes prediction on the points belonging to this candidate track, for example on the basis of the candidate track found already and according to the presently-used track model to predict the other hit points by interpolation or extrapolation. If extra points have been found, then they are added to the candidate track, otherwise, after trying for a certain times (the number of times depends on the tolerance of detector failure allowed by the algorithm), this candidate track is discarded. In order to find a candidate track, the localized method always makes the resultless trial, thus the different combinations may use the same point, as the number of points increases, the growth rate of computing time will be far greater than the rate of linear growth. If all objects (hit points) appear in the algorithm after the same manner, then it is called the overall method. This algorithm produces a track list, by which the track can be more easily found from the original data, this algorithm can be considered as a kind of general transformation of the whole data set of event coordinates or space points. In principle, the computing time of the overall method will be directly proportional to the number of hit points in an event. For this particular detector, when a high-energy electron enters into the silicon microstrip detector, it will lose a little energy because of the ionization loss, but in the STK there exist tungsten plates of totally 3 mm thickness, a very large part of electrons will produce electromagnetic showers[11,12] in the tungsten plates, and produce even more electrons which will be detected after being converted into electric signals. Many clusters may be produced in the several succeeding silicon microstrip layers, of course, the interaction of the back-scattered particle with the upper silicon microstrip layer may also produce a cluster, but it is much smaller compared to the electromagnetic shower. We select the position of cluster centroid as a point in the process of data treatment. When a high-energy proton enters into the silicon microstrip detector, it happens mainly the ionization loss, although the thickness of tungsten plate is small, but the density is rather high, a small number of hadron showers may happen as well, but the probability is quite small. According to the physical principle mentioned above, combining with the restrictive conditions, the ergodic search is performed for all possible tracks. In which one method is: in the ergodic search, we select arbitrarily two points in the restricted range to fit a straight
462
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
line, and judge whether there are additional points in the other layers close to this straight line (as shown in Fig.6, if the distance from a point to the straight line ls is less than 1 mm, then this point is considered to be close to the straight line), so that to find out that straight line which has most points close to it with the minimum variance σ 2 , then using the points on this straight line and the points close to this straight line, we fit again a straight line, this newly fitted straight line is just the reconstructed track in the silicon microstrip detector that we are looking for. In order to reduce the amount of calculation, and to improve the accuracy, the finally selected candidate track contains at least 3 points, ideally 6 points, in this case even if there is a bad point, there are still 5 points to participate the fitting, if the selected points are less than 3, even though 2 points can make a straight line, but the appearance probability of this case is very small[13] , almost negligible. l2
lc
l1
ls
5°
Fig. 6 The geometric range for searching the cluster
In the searching process the following two aspects are mainly taken into consideration: one is that the track should be as long as possible, and should be a straight line as possible, speaking popularly, the selected points are as many as possible, and the variance is as small as possible. The second is that it should match most the direction of the BGO track, or be consistent as possible, namely both the slope and intercept are as approximate as possible, but in a practical treatment the comparison on the slope and intercept is rather difficult. Hence, we adopt a more reasonable and convenient method, as shown in Fig.7, in this figure the two straight lines express respectively the track fitted and reconstructed by using the BGO information, and a possible track obtained by using the silicon microstrip information and after searching all possible points, their points of intersection with the silicon microstrip layers are respectively X1 , X2 , · · · , X6 and X1 , X2 , · · · , X6 . We can obtain very easy that when n (Xi − Xi )2 (3) i=1
is smaller, the two straight lines are closer to each other. Namely, when it is minimum, the two straight lines match best. Another method is: as arbitrary two points can determine a straight line, we select arbitrarily two points in the restricted range to fit a straight line,
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
463
and the straight line which deviates the track reconstructed in the BGO calorimeter by more than 5◦ is deleted, because basically most points (possible track) are included within a range of 5◦ . Then, to judge whether there are points on other layers which are close to the straight line, and find out the straight line with most points close to it and with the minimum χ2 -value, using the points on this straight line and the points close to it, we make again the linear fitting, the newly fitted straight line is just the reconstructed track in the silicon microstrip detector that we are looking for. The track in BGO Silicon Calorimeter microstrip X1 X1′ X 2
X3′ X4′ X5′ X6′
X2′
X3 X4
X5
X6
Z
X
The candidate track in the STK Fig. 7 The projection of the candidate track in the STK on the XZ plane
4.3
Track Fitting
4.3.1 Least squares linear fitting with energy weights For the known data points (xi , yi ) (i = 1, 2, · · · , n), to make the energy weighting xi = Exi xi , yi = Eyi yi , and fitting with y = b + a x, we obtain a and b as follows n n n n 2 i=1 yi i=1 xi − i=1 xi i=1 xi yi b = , n n 2 2 n i=1 xi − ( i=1 xi ) (4) n n n n i=1 xi yi − i=1 xi i=1 yi n , a = n n i=1 x2i − ( i=1 xi )2 in which Exi and Eyi express respectively the energy values corresponding to the xi -th and yi -th clusters or crystals, a and b are respectively the slope and intercept of the straight line. The results of reconstruction are shown as Fig.8 and Fig.9. Fig.8 shows the result of track reconstruction for an electron of 100 GeV incident from the direction of 30◦ , in which the solid line is the fitting result of the BGO centroid method, the dashed line is the restriction given by the BGO calorimeter, and the dot-dashed line is the best-fit result of the track in the STK, which is the projection on the XZ plane, the abscissa expresses the centroid position of hit points along the Z axis, and the ordinate indicates the centroid position of hit points along the X axis. Fig.9 shows the result of track reconstruction for a proton of 400 GeV with the incident angle of 0◦ .
464
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
Fig.10 verifies very well the effectiveness of track reconstruction, the upper left panel expresses the distribution of the angle θx between the projection of the particle track in the BGO calorimeter on the XZ plane and the Z axis with the event accumulation; the lower right panel shows the distribution of the angle θy between the projection of the particle track in the BGO calorimeter on the Y Z plane and the Z axis with the event accumulation; in the lower left panel, the abscissa expresses the angle θx between the projection of particle track on the XZ plane and the Z axis, and the ordinate expresses the angle θy between the projection of particle track on the Y Z plane and the Z axis. The upper right panel indicates the distribution of the angle θ between the particle track and the Z axis with the event accumulation. In the figure, for the ideal result of track reconstruction of a vertically incident particle, the central position of the angle distribution should be about 0◦ , but for the reconstructed result of a vertically incident proton of 400 GeV, as shown in Fig.10, in the Y direction it is basically positioned at the zero point, but in the X direction it deviates about 0.35◦ from the zero point, this may be caused by two reasons: (1) because of the installation error of the detector and some systematic error of the measurement system, the particle is not exactly incident from the direction of 0◦ ; (2) our algorithm remains to be improved, it is failed to find the real track. From the upper left and lower right panels we can find that there are tails deviated from the central position, these are basically the false tracks occurred in the reconstruction, but not completely deleted by us. And from the lower left panel we can see also some scattered points deviated from the central position. 300 250 200 150 100 50 0 -50 -100 -200
-100
0
100
200
300
400
500
Fig. 8 The trajectory of an electron with a momentum of 100 GeV and an incident angle of 30◦ reconstructed by combining the BGO calorimeter and STK
4.3.2 The linear fitting based on Kalman filtering Kalman filtering is a kind of linear filtering and prediction method proposed first by the Hungarian mathematician Rudof Email Kalman in 1960[14] . In the linear fitting, the least squares method is often used[15,16] , for the linear fitting with pseudo-points included, be-
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
465
cause a part of pseudo-points have a large deviation from the fitted straight line, the fitted track may deviate from the real track, some pseudo-points are caused artificially, in this case we may delete these pseudo-points before fitting. But to delete the pseudo-points by experience or intuition is not believable enough, they should be judged by some effective methods. Kalman filtering is a kind of statistical estimation method, because of its property of minimum unbiased variance, it can remove the random errors in the experiment to obtain the information more close to the reality, it has been widely applied to different fields sofar. For the solution of a major part of problems, it is a kind of most useful algorithm with the highest efficiency.
300
200
100
0
-100
-200
-100
0
100
200
300
400
500
Fig. 9 The trajectory of a proton with a momentum of 400 GeV and an incident angle of 0◦ reconstructed by combining the BGO calorimeter and STK
Figs.11∼12 are the scatter plots of the slopes and intercepts of the tracks after a cosmic muon hits the detector as derived respectively by using the energy-weighted least squares method and Kalman filtering. In which, the abscissa of Fig.11 expresses the slope (inclination) of the track obtained by the energy-weighted least squares linear fitting, in units of degree; the ordinate expresses the slope (inclination) of the track obtained by the Kalman filtering method, in units of degree; the abscissa of Fig.12 expresses the intercept of the track obtained by the energy-weighted least squares linear fitting, in units of mm; the ordinate expresses the intercept of the track obtained by the Kalman filtering method, in units of mm. From Figs.11∼12 we can see clearly that both distributions of slopes and intercepts are rather concentrated, almost all are distributed along a straight line with the slope equal to 1, indicating the equivalent reconstruction effectiveness of the Kalman filtering and energy-weighted least squares fitting, and that the difference between the two is not very large. From the data produced when a 100 GeV electron or a 400 GeV proton is vertically
466
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
10000
10000
8000
8000 Entries
Entries
incident to the detector, we randomly select two events, the parameters of the two tracks fitted by using the energy-weighted least squares method and Kalman filtering method are compared in Table 1, respectively.
6000 4000
6000 4000
2000
2000
5 4 3 2 1 0 1 2 3 4 5 4 3 2 1 0
1 2 θx /(°)
3
4
5
6
0
1
2
3
4
5 6 θ /(°)
7
8
9 10
0 1 θy /(°)
2
3
4
14000 12000 Entries
θy /(°)
0 4 3 2 1 0
1 2 θx /(°)
3
4
5
10000 8000 6000 4000 2000 0 5 4 3 2 1
6
5
Fig. 10 The distribution of track directions in the STK
60 40 20 0 -20 -40 -60 -60
-40
-20
0
20
40
60
Fig. 11 The scatter plot of the slopes (inclinations) fitted with two methods
The results of track reconstruction are given in Fig.13 and Fig.14, in which the solid line is the track reconstructed by using the Kalman filtering method, the dotted line is the result of energy-weighted least squares linear fitting, the dashed line is the result reconstructed by using the shower in the BGO calorimeter, and the black spots indicate the hit points in
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
467
the detector. We have made the statistics on the results of track reconstruction at multiple energy bands, and find that in the low-energy region, the fitting result based on the Kalman filtering is obviously superior to that of the least squares fitting, as shown in Fig.13, the solid line and the dotted line are not coincident, and the solid line is more approximate to the real track. In the high-energy region, the results of the energy-weighted least squares method and Kalman filtering method are not very different, as shown in Fig.14, the solid line and dotted line are basically coincident, and after the position calibration of the detector, they are all approximate to the real track. 500 400 300 200 100 0 -100 -200 -300 -400 -500 -500
-400 -300 -200 -100
0
100
200
300
400
500
Fig. 12 The scatter plot of the intercepts fitted with two methods
Table 1 The parameters of track fitting 100 GeV electron
Parameter
400 GeV proton
Slope
Intercept/mm
Slope
Intercept/mm
Least square method
0.01974
128.336
0.00034237
178.704
Kalman filter
0.01078
129.237
0.0010725
178.759
Method
5.
ERROR ANALYSIS
There are many factors to cause the error in the result of track reconstruction, the main factors can be described as follows. (1) In the beam experiment, the detector position itself is approximate, not very accurate, or some parts of the detector may deviate from the central position, and the deviated directions may be different, this will influence on the selection of coordinate system, and cause the error of position calculation of a certain part of the detector relative to the coordinate axis or coordinate origin. It requires to make the position calibration before the data processing, ensuring that the calculations are performed in the same coordinate system.
468
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
(2) Because of the material of the detector, when a particle hits the detector, many secondary particles may be produced, including some charged particles, we can not judge that the cluster is composed of what kinds of particles. Of course, it is possible that a cluster may be produced also by the interaction of a recoil particle with the detector. (3) In the data processing, the centroid position of the cluster is selected as the hit point, but the position that a particle actually hits the detector is not certainly the centroid position, there is a deviation existed. 250 200 150 100 50 0 -200
-100
0
100
200
300
400
500
Fig. 13 The reconstructed trajectory of an electron with a momentum of 100 GeV and an incident angle of 0◦
300 200 100 0 -100 -200 -300 -200
-100
0
100
200
300
400
500
Fig. 14 The reconstructed trajectory of a proton with a momentum of 400 GeV and an incident angle of 0◦
6.
SUMMARY
The track reconstruction has an important significance for the detection of dark matter particles, to recover the tracks of incident particles will help us to obtain the information
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
469
about the type, energy, intensity, and incident direction of incident particles, and to analyze further the generation of particles and the direction of the source. If these particles are produced by the decay or annihilation of dark matter, we can learn about the distribution of dark matter, as well as the correlation with the distribution of other kind matter. The most important points for the track reconstruction are: the deletion of pseudo-points, track search, track fitting, and the relevant corrections. Several key points should be taken into consideration in the process of track reconstruction, for example the detector performance (mainly the geometric configuration and accuracy (resolution)), energy loss, multiple scattering, etc., so that we can build a mathematic model to describe the particle track with a fairly good accuracy. In order to verify the accuracies of the model and algorithm, before analyzing the actual experimental data, the whole reconstruction process should be tested by using the simulated data, and the real values of estimated parameters in the simulated data should be known already, in favor of the test on the normalized residuals. The exact errors of the fitted parameters should be investigated carefully, which may be influenced by the pseudo-points or pseudo-tracks, as well as the method of their deletion, and the position error of the detector installation. This paper adopts mainly the energy-weighted least squares method and the straight line fitting method based on the Kalman filtering to reconstruct the incident particle track, by comparisons, it is found that the result (error) of track reconstruction is related with the energy of the incident particle. Generally, the higher the incident energy, the smaller the error of the result of track reconstruction. Further more, we find that in the high-energy region, the result obtained by the energy-weighted least squares linear fitting and the result obtained after the prediction, filtering and smoothing of the Kalman filter are basically equivalent, but in the low-energy region, the Kalman filtering has a relatively apparent superiority. The Kalman filter is a kind of powerful, flexible, and effective tool, it can be used for not only the track fitting, but also the calculations of optimum extrapolation and interpolation, the search and deletion of pseudo-points, the error adjustment and detector calibration using the actual data, as well as the merge of track segments. However, the Kalman filter has also its own defects, some times it is difficult to find a suitable initial value of Kalman filtering, if the initial angular resolution and momentum resolution of the detector are very poor, and the next detector is not very near, the first-step prediction may have the following problems: (1) It is very difficult to calculate precisely the amount of matter to be penetrated by the real particle track, because the predicted trajectory may be completely different from the real trajectory. (2) If the initial value deviates far from the real track, the linearity approximation of the track model will be invalid. (3) When the predicted track does not intersect with the surface of the next detector, the prediction procedure will be unable to perform.
470
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
A method to solve this kind of problems is to take the reference track derived by the pattern recognition method as the initial point of linear approximation. Namely, to calculate first the intersections of the reference track with all surfaces (the so-called reference state), then to make prediction, filtering, and smoothing by using the difference between the state and the reference state. Besides, the position measurement during the instrument installation is not perfectly accurate, the calibration of the instrument will have a certain influence on the result of track reconstruction, hence, in order to improve the effectiveness of the track reconstruction, it is also necessary to make calibration on the instrument. References 1
Chang J. ChJSS, 2014, 34, 550
2
Dong Y. F., Study of Data Acquistion System for the Silicon Detector on the Dark Matter Detection Satellite, Bejing: Institute of High-energy Physics of Chinese Academy of Sciences, 2015
3
Zhang Y. L., Study of Electromagnetic Calorimeter for Dark Matter Detection, Hefei: The University of Scence and Technology of China, 2011
4
Xie M. G., Guo J. H., Wu J., et al., AcASn, 2014, 55, 170
5
Xie M. G., Guo J. H., Wu J., et al., ChA&A, 2015, 39, 129
6
Guo J. H., Cai M. S., Hu Y. M., et al., AcASn, 2012, 53, 72
7
Zhang L., Guo J. H., Zhang Y. Q., AcASn, 2014, 55, 522
8
Zhang L., Guo J. H., Zhang Y. Q., ChA&A, 2015, 39, 380
9
Fruhwirth R., Regle M., High-energy Physics Data Analysis, Translatged by Zhu Y. S., Hefei: Press of the University of Science and Technology of China, 2011
10 11
Liu J., Mao Z. P., Li W. G., et al., High Energy Physics and Nuclear Physics, 1998, 22, 587 Wang X. L., Particle Detection Technique, Hefei: The University of Science and Technology of China, 2009
12
Mao Z. P,, Wang T. J., Wang S. Q., HEPNP, 1993, 17, 193
13
Grupen C., Shwartz B., Spieler H., Particle Detectors, Cambridge: Cambridge University Press, 2011
14
Kalman R. E., Transactions of the ASME-Journal of Basic Engineering, 1960, 35
15
Wang J. R., Instrument and Meter for Automation, 2013, 140
16
Gao F., Dong H. Q., System Engineering and Electronic Technology, 2006, 28, 775
Study of Track Reconstruction for DAMPE† LU Tong-suo1,2,3
LEI Shi-jun1,2
CHANG Jin1,2 1 2
ZANG Jing-jing1,2
WU Jian1,2
Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008
Key Laboratory of Dark Matter and Space Astronomy, Chinese Academy of Sciences, Nanjing 210008 3
University of Chinese Academy of Sciences, Beijing 100049
Abstract The Dark Matter Particle Explorer (DAMPE) is aimed to study the existence and distribution of dark matter via the observation of high-energy particles in space with a large energy bandwidth, high energy resolution, and high spatial resolution. The track reconstruction is to restore the positions and angles of the incident particles using the multiple observations of different channels at different positions, and its accuracy determines the angular resolution of the explorer. The track reconstruction is mainly based on the observations of two sub-detectors, namely, the Silicon Tracker (STK) and the BGO (Bi4 Ge3 O12 ) calorimeter. In accordance with the design and structure of the two sub-detectors, and using the data collected during the beam tests and ground tests of cosmic rays, we discuss in detail the method of track reconstruction for the DAMPE, which includes mainly three basic procedures: the selection of track hits, the fitting of track hits, and the judgement of the optimal track. Since an energetic particle most probably leaves multiple hits in different layers of the STK and BGO crystals, we first provide a method to obtain a rough track in the BGO calorimeter by the centroid method, and hereby to constrain the track hits in the STK. Then for the selected one group of possible track hits in the STK, we apply two different algorithms, the Kalman filter and the least square linear fitting, to fit these track hits. The consistency of the results obtained independently via the two algorithms confirms the validity of our track reconstruction results. Finally, several criteria for picking out the most possible track among all the tracks found in the reconstruction by combining the results †
Supported by National Natural Science Foundation (11303105, 11303106, 11303107)
A translation of Acta Astron. Sin. Vol. 57, No. 3, pp. 353–365, 2016
Received 2015–10–16; revised version 2015–11–16
[email protected]
[email protected]
0275-1062/17/$-see front matter © 2017 Elsevier B.V. All rights reserved. doi:10.1016/j.chinastron.2017.08.012
456
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
of the BGO calorimeter and STK are discussed. Using the track reconstruction method proposed in this article and the beam test data, we confirm that the angular resolution of the DAMPE satisfies its design requirement. Key words instrumentation: detectors—planets and satellites: detection— gamma rays: ISM— cosmology: dark matter—methods: data analysis 1.
INTRODUCTION
Dark matter is a puzzle in astronomy and the development of modern astrophysics, now the detection of dark matter is during the period of rapid development. In China, aiming at the dark matter, the Dark Matter Particle Explorer (DAMPE) was launched into space on 17th December 2015, it adopts a 3-dimensional total-absorption calorimeter of large dynamic range, and has a series of advantages, such as the large energy range, high energy resolution, and the strong ability of back ground suppression, etc. According to the requirement of its physical design and the result of simulations, the DAMPE is composed of 4 sub-detectors, from top to bottom, they are the Plastic Scintillator Detector (PSD), Silicon Tracker (STK), BGO calorimeter, and Neutron Detector (NUD), as shown by Fig.1. The effective detection area of the PSD subsystem is 820 mm×820 mm, the whole subsystem consists of two layers of totally 82 plastic scintillator modules perpendicular to each other along the X and Y directions, in which the size of 78 plastic scintillator modules is 884 mm×28 mm×10 mm, and the size of other four modules is 884 mm×25×10 mm. On the both ends of the PSD unit, photomultipliers are adopted to convert the optical signals into electric signals for further treatment. The STK has 6 major layers, each major layer consists respectively of two layers of silicon microstrip detectors perpendicular to each other along the X and Y directions, mainly for measuring the directions and tracks of incident particles. A piece of tungsten plate of 1 mm thickness is placed respectively upon the second, third, and fourth major layers, to convert Gamma photons into positive and negative electrons. Considered that after the photon conversion, it is necessary to measure the tracks of the produced secondary particles, in order to distinguish the charged particles, especially distinguish the Gamma photons and electrons[1] , under the tungsten plate there must be multiple layers of silicon microstrip detectors[2] . AS shown by Fig.2, the main body of the BGO calorimeter is made of 308 pieces of BGO crystals, its detection area is 60 cm×60 cm. The BGO calorimeter is composed of totally 7 major layers, each major layer consists of the X and Y two minor layers. And each minor layer has 22 detector units composed of 22 pieces of crystals. The size of each detector unit (BGO crystal) is 2.5 cm×2.5 cm×60 cm. On the two ends of the detector unit, a R5610A-01 photomultiplier (PMT) is adopted to convert the optical signals into electric signals[3−6] . The NUD is positioned on the bottom of the DAMPE, it adopts the B-doped plastic
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
457
scintillator detector with a thickness of 1 cm (BC454, made by the Saint Gobian Company). The detector is composed of the plastic scintillator with a profile size of 693 mm×693 mm, which is divided into 4 independent squares, and one angle of each square is cut off for coupling with one PMT readout to form a complete detection plane. It is designed for measuring the secondary particles produced by the interactions of the hadrons in cosmic rays (mainly protons) with the matter on the upper layer of the NUD, according to the energy precipitation of these protons in the detector, the types of incident particles can be judged, and in cooperation with the BGO calorimeter, the protons and electrons can be further distinguished[7,8] .
Fig. 1 A sketch of 4 sub-detectors
2.
Fig. 2 The schematic diagram of the BGO calorimeter
PROCEDURES OF TRACK RECONSTRUCTION
Considered the geometric configuration and accuracy of the detector, in order to obtain a mathematic model to describe approximately the particles, we adopt the track fitting method shown in Fig.3. 3.
TRACK RECONSTRUCTION OF PARTICLES IN THE BGO CALORIMETER
The track reconstruction of particles in the BGO calorimeter is based on the following consideration: when the incident particles pass through the calorimeter, because of the interactions (mainly the electromagnetic interaction and strong interaction) of particles with the medium atoms, the particles will deposit energies in the calorimeter, according to the shower shape, we can judge approximately the track direction. The practical realization is as follows: at first, to select a BGO layer, by setting an energy threshold to get rid of the background noise, according to multiple statistical results, we preliminarily set the threshold
458
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
as 2%, only when the energy of a BGO layer exceeds 2% total energy, can this layer be taken as the reference layer to join in the track reconstruction in the BGO calorimeter. Secondly, to select the BGO crystals in each BGO layer to participate in the track reconstruction, we select the crystal with the greatest energy deposition and the two crystals in its left and right sides, totally 3 crystals, to calculate the centroid of energy deposition of this layer, the details are shown as Eq.(1). If in a layer the crystal with the greatest energy deposition is positioned on the edge, namely there is no crystal in its left or right side, in order to reduce errors, we select only this crystal with the greatest energy deposition to join in the track reconstruction; if in a layer, the crystal with the greatest energy deposition is not positioned on the edge, but the energies of the crystals in its right and left sides have a large difference, namely the energy of one side is negligible in respect to the other side, we still take the two crystals with the higher energies to participate the reconstruction, because the situation that a particle penetrates through simultaneously two neighboring crystals and gets out from the underneath is possible. Then, to make linear fitting: the energy-weighted linear fitting, the energy used for weighting is the energy summation of all the crystals participating in the reconstruction. Emax−1 xmax−1 + Emax xmax + Emax+1 xmax+1 = xc , Emax−1 + Emax + Emax+1
(1)
In Eq.(1), xc , xmax , and Emax indicate respectively the position of energy centroid of a layer in the BGO calorimeter, the position and energy of the crystal with the greatest energy in this layer. xmax−1 , xmax+1 , Emax−1 , and Emax+1 express respectively the positions and energies of the crystals in the left and right of xmax . Fig.4 shows the result of track reconstruction of a vertically incident proton of 400 GeV. This figure is the projective diagram on the XZ plane, the several lines of dots in the right side are the energy centroid positions of the crystals actually being hit in the calorimeter, the length of each transverse error bar indicates the magnitude of energy deposition, the straight line is the fitted particle track. In this figure the 6 points (X <0) belong to the cluster in the STK. 4. 4.1
TRACK RECONSTRUCTION OF PARTICLES IN THE STK
The Selection of Cluster
As the size of the BGO crystal (2.5 cm×2.5 cm×60 cm) is relatively large, so some errors exist unavoidably in the process of track reconstruction, this is caused by the system itself, we can use this error range to constrain the selection of cluster in the STK. Fig.5 shows the distribution of particle track direction reconstructed by using the BGO calorimeter information, which can be taken as a reference for the particle track reconstruction in the STK. The upper left panel displays the distribution of the angle θx between the projection of the particle track in the BGO calorimeter on the XZ plane and the Z axis,
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
459
Set the thresholdˈfind out the position of energy center
Fitting the track in BGO calorimeter
Selecting the clusters
Selecting the candidate tracks
Find out the best track
Kalman filter
The track Fig. 3 Flow chart for reconstructing the incident particle trajectory
300
X/mm
200 100 0 100 200 200
100
0
100 Z/mm
200
300
400
Fig. 4 The result of track reconstruction of a proton vertically incident to the top surface of the DAMPE with a momentum of 400 GeV
460
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
as a function of event accumulation; the lower right panel expresses the distribution of the angle θy between the projection of the particle track in the BGO calorimeter on the Y Z plane and the Z axis, as a function of event accumulation; in the lower left panel the abscissa expresses the angle θx between the projection of the particle track in the BGO calorimeter on the XZ plane and the Z axis, the ordinate expresses the angle θy between the projection of the particle track in the BGO calorimeter on the Y Z plane and the Z axis. The upper right panel shows the distribution of the angle θ between the particle track and the Z axis as a function of event accumulation. According to the calculation, in order to reduce the errors, the deviation of the slope (inclination) of the reconstructed track in the STK from that in the BGO calorimeter should be within ±5◦ , from Fig.5 we can find that most hit points are distributed in a 5◦ range, and the intercept is taken to be ±25 mm, namely the search of particle track in the STK is restricted in this range. We assume that the equation of projection of the reconstructed track in the BGO calorimeter on the XZ plane is l c : z = ax + b, as shown by Fig.6. Then, the equations of the straight lines which restrict the cluster search range are: l2 : z = tan(arctan a + 5)x + 25 + b ,
(2)
l1 : z = tan(arctan a − 5)x − 25 + b ,
5 4 3 2 1 0 1 2 3 4 5 4 3 2 1
1000 Entries
800 600 400 200 0
1 2 θx /(°)
3
4
5
6
0
Entries
1600 1400 1200 1000 800 600 400 200 0 4 3 2 1
θy /(°)
Entries
in which, a and b are respectively the slope (inclination) and intercept of the equation of projection.
0
1 2 θx /(°)
3
4
5
6
2
4
6
8
1800 1600 1400 1200 1000 800 600 400 200 0 5 4 3 2 1
Fig. 5 The track reconstruction in the BGO calorimeter
10 12 14 16 18 20 θ /(°)
0 1 θy /(°)
2
3
4
5
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
4.2
461
Track Search
Track search is a kind of pattern recognition, belonging to the scope of artificial intelligence[9,10] , its task is to divide a group of experimentally measured points into two types: the first type can be divided into many subtypes, the points in each subtype are all caused by one and the same particle, namely, corresponding to one track; the second type indicates those points which can not be attributed to any particle under sufficient confidence level, namely the noise. At present, the method of track recognition is of two kinds, namely the localized method and overall method. In the localized method, each time only one candidate track is selected, typically it starts from a few points (the initialization of candidate track), then it makes prediction on the points belonging to this candidate track, for example on the basis of the candidate track found already and according to the presently-used track model to predict the other hit points by interpolation or extrapolation. If extra points have been found, then they are added to the candidate track, otherwise, after trying for a certain times (the number of times depends on the tolerance of detector failure allowed by the algorithm), this candidate track is discarded. In order to find a candidate track, the localized method always makes the resultless trial, thus the different combinations may use the same point, as the number of points increases, the growth rate of computing time will be far greater than the rate of linear growth. If all objects (hit points) appear in the algorithm after the same manner, then it is called the overall method. This algorithm produces a track list, by which the track can be more easily found from the original data, this algorithm can be considered as a kind of general transformation of the whole data set of event coordinates or space points. In principle, the computing time of the overall method will be directly proportional to the number of hit points in an event. For this particular detector, when a high-energy electron enters into the silicon microstrip detector, it will lose a little energy because of the ionization loss, but in the STK there exist tungsten plates of totally 3 mm thickness, a very large part of electrons will produce electromagnetic showers[11,12] in the tungsten plates, and produce even more electrons which will be detected after being converted into electric signals. Many clusters may be produced in the several succeeding silicon microstrip layers, of course, the interaction of the back-scattered particle with the upper silicon microstrip layer may also produce a cluster, but it is much smaller compared to the electromagnetic shower. We select the position of cluster centroid as a point in the process of data treatment. When a high-energy proton enters into the silicon microstrip detector, it happens mainly the ionization loss, although the thickness of tungsten plate is small, but the density is rather high, a small number of hadron showers may happen as well, but the probability is quite small. According to the physical principle mentioned above, combining with the restrictive conditions, the ergodic search is performed for all possible tracks. In which one method is: in the ergodic search, we select arbitrarily two points in the restricted range to fit a straight
462
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
line, and judge whether there are additional points in the other layers close to this straight line (as shown in Fig.6, if the distance from a point to the straight line ls is less than 1 mm, then this point is considered to be close to the straight line), so that to find out that straight line which has most points close to it with the minimum variance σ 2 , then using the points on this straight line and the points close to this straight line, we fit again a straight line, this newly fitted straight line is just the reconstructed track in the silicon microstrip detector that we are looking for. In order to reduce the amount of calculation, and to improve the accuracy, the finally selected candidate track contains at least 3 points, ideally 6 points, in this case even if there is a bad point, there are still 5 points to participate the fitting, if the selected points are less than 3, even though 2 points can make a straight line, but the appearance probability of this case is very small[13] , almost negligible. l2
lc
l1
ls
5°
Fig. 6 The geometric range for searching the cluster
In the searching process the following two aspects are mainly taken into consideration: one is that the track should be as long as possible, and should be a straight line as possible, speaking popularly, the selected points are as many as possible, and the variance is as small as possible. The second is that it should match most the direction of the BGO track, or be consistent as possible, namely both the slope and intercept are as approximate as possible, but in a practical treatment the comparison on the slope and intercept is rather difficult. Hence, we adopt a more reasonable and convenient method, as shown in Fig.7, in this figure the two straight lines express respectively the track fitted and reconstructed by using the BGO information, and a possible track obtained by using the silicon microstrip information and after searching all possible points, their points of intersection with the silicon microstrip layers are respectively X1 , X2 , · · · , X6 and X1 , X2 , · · · , X6 . We can obtain very easy that when n (Xi − Xi )2 (3) i=1
is smaller, the two straight lines are closer to each other. Namely, when it is minimum, the two straight lines match best. Another method is: as arbitrary two points can determine a straight line, we select arbitrarily two points in the restricted range to fit a straight line,
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
463
and the straight line which deviates the track reconstructed in the BGO calorimeter by more than 5◦ is deleted, because basically most points (possible track) are included within a range of 5◦ . Then, to judge whether there are points on other layers which are close to the straight line, and find out the straight line with most points close to it and with the minimum χ2 -value, using the points on this straight line and the points close to it, we make again the linear fitting, the newly fitted straight line is just the reconstructed track in the silicon microstrip detector that we are looking for. The track in BGO Silicon Calorimeter microstrip X1 X1′ X 2
X3′ X4′ X5′ X6′
X2′
X3 X4
X5
X6
Z
X
The candidate track in the STK Fig. 7 The projection of the candidate track in the STK on the XZ plane
4.3
Track Fitting
4.3.1 Least squares linear fitting with energy weights For the known data points (xi , yi ) (i = 1, 2, · · · , n), to make the energy weighting xi = Exi xi , yi = Eyi yi , and fitting with y = b + a x, we obtain a and b as follows n n n n 2 i=1 yi i=1 xi − i=1 xi i=1 xi yi b = , n n 2 2 n i=1 xi − ( i=1 xi ) (4) n n n n i=1 xi yi − i=1 xi i=1 yi n , a = n n i=1 x2i − ( i=1 xi )2 in which Exi and Eyi express respectively the energy values corresponding to the xi -th and yi -th clusters or crystals, a and b are respectively the slope and intercept of the straight line. The results of reconstruction are shown as Fig.8 and Fig.9. Fig.8 shows the result of track reconstruction for an electron of 100 GeV incident from the direction of 30◦ , in which the solid line is the fitting result of the BGO centroid method, the dashed line is the restriction given by the BGO calorimeter, and the dot-dashed line is the best-fit result of the track in the STK, which is the projection on the XZ plane, the abscissa expresses the centroid position of hit points along the Z axis, and the ordinate indicates the centroid position of hit points along the X axis. Fig.9 shows the result of track reconstruction for a proton of 400 GeV with the incident angle of 0◦ .
464
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
Fig.10 verifies very well the effectiveness of track reconstruction, the upper left panel expresses the distribution of the angle θx between the projection of the particle track in the BGO calorimeter on the XZ plane and the Z axis with the event accumulation; the lower right panel shows the distribution of the angle θy between the projection of the particle track in the BGO calorimeter on the Y Z plane and the Z axis with the event accumulation; in the lower left panel, the abscissa expresses the angle θx between the projection of particle track on the XZ plane and the Z axis, and the ordinate expresses the angle θy between the projection of particle track on the Y Z plane and the Z axis. The upper right panel indicates the distribution of the angle θ between the particle track and the Z axis with the event accumulation. In the figure, for the ideal result of track reconstruction of a vertically incident particle, the central position of the angle distribution should be about 0◦ , but for the reconstructed result of a vertically incident proton of 400 GeV, as shown in Fig.10, in the Y direction it is basically positioned at the zero point, but in the X direction it deviates about 0.35◦ from the zero point, this may be caused by two reasons: (1) because of the installation error of the detector and some systematic error of the measurement system, the particle is not exactly incident from the direction of 0◦ ; (2) our algorithm remains to be improved, it is failed to find the real track. From the upper left and lower right panels we can find that there are tails deviated from the central position, these are basically the false tracks occurred in the reconstruction, but not completely deleted by us. And from the lower left panel we can see also some scattered points deviated from the central position. 300 250 200 150 100 50 0 -50 -100 -200
-100
0
100
200
300
400
500
Fig. 8 The trajectory of an electron with a momentum of 100 GeV and an incident angle of 30◦ reconstructed by combining the BGO calorimeter and STK
4.3.2 The linear fitting based on Kalman filtering Kalman filtering is a kind of linear filtering and prediction method proposed first by the Hungarian mathematician Rudof Email Kalman in 1960[14] . In the linear fitting, the least squares method is often used[15,16] , for the linear fitting with pseudo-points included, be-
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
465
cause a part of pseudo-points have a large deviation from the fitted straight line, the fitted track may deviate from the real track, some pseudo-points are caused artificially, in this case we may delete these pseudo-points before fitting. But to delete the pseudo-points by experience or intuition is not believable enough, they should be judged by some effective methods. Kalman filtering is a kind of statistical estimation method, because of its property of minimum unbiased variance, it can remove the random errors in the experiment to obtain the information more close to the reality, it has been widely applied to different fields sofar. For the solution of a major part of problems, it is a kind of most useful algorithm with the highest efficiency.
300
200
100
0
-100
-200
-100
0
100
200
300
400
500
Fig. 9 The trajectory of a proton with a momentum of 400 GeV and an incident angle of 0◦ reconstructed by combining the BGO calorimeter and STK
Figs.11∼12 are the scatter plots of the slopes and intercepts of the tracks after a cosmic muon hits the detector as derived respectively by using the energy-weighted least squares method and Kalman filtering. In which, the abscissa of Fig.11 expresses the slope (inclination) of the track obtained by the energy-weighted least squares linear fitting, in units of degree; the ordinate expresses the slope (inclination) of the track obtained by the Kalman filtering method, in units of degree; the abscissa of Fig.12 expresses the intercept of the track obtained by the energy-weighted least squares linear fitting, in units of mm; the ordinate expresses the intercept of the track obtained by the Kalman filtering method, in units of mm. From Figs.11∼12 we can see clearly that both distributions of slopes and intercepts are rather concentrated, almost all are distributed along a straight line with the slope equal to 1, indicating the equivalent reconstruction effectiveness of the Kalman filtering and energy-weighted least squares fitting, and that the difference between the two is not very large. From the data produced when a 100 GeV electron or a 400 GeV proton is vertically
466
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
10000
10000
8000
8000 Entries
Entries
incident to the detector, we randomly select two events, the parameters of the two tracks fitted by using the energy-weighted least squares method and Kalman filtering method are compared in Table 1, respectively.
6000 4000
6000 4000
2000
2000
5 4 3 2 1 0 1 2 3 4 5 4 3 2 1 0
1 2 θx /(°)
3
4
5
6
0
1
2
3
4
5 6 θ /(°)
7
8
9 10
0 1 θy /(°)
2
3
4
14000 12000 Entries
θy /(°)
0 4 3 2 1 0
1 2 θx /(°)
3
4
5
10000 8000 6000 4000 2000 0 5 4 3 2 1
6
5
Fig. 10 The distribution of track directions in the STK
60 40 20 0 -20 -40 -60 -60
-40
-20
0
20
40
60
Fig. 11 The scatter plot of the slopes (inclinations) fitted with two methods
The results of track reconstruction are given in Fig.13 and Fig.14, in which the solid line is the track reconstructed by using the Kalman filtering method, the dotted line is the result of energy-weighted least squares linear fitting, the dashed line is the result reconstructed by using the shower in the BGO calorimeter, and the black spots indicate the hit points in
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
467
the detector. We have made the statistics on the results of track reconstruction at multiple energy bands, and find that in the low-energy region, the fitting result based on the Kalman filtering is obviously superior to that of the least squares fitting, as shown in Fig.13, the solid line and the dotted line are not coincident, and the solid line is more approximate to the real track. In the high-energy region, the results of the energy-weighted least squares method and Kalman filtering method are not very different, as shown in Fig.14, the solid line and dotted line are basically coincident, and after the position calibration of the detector, they are all approximate to the real track. 500 400 300 200 100 0 -100 -200 -300 -400 -500 -500
-400 -300 -200 -100
0
100
200
300
400
500
Fig. 12 The scatter plot of the intercepts fitted with two methods
Table 1 The parameters of track fitting 100 GeV electron
Parameter
400 GeV proton
Slope
Intercept/mm
Slope
Intercept/mm
Least square method
0.01974
128.336
0.00034237
178.704
Kalman filter
0.01078
129.237
0.0010725
178.759
Method
5.
ERROR ANALYSIS
There are many factors to cause the error in the result of track reconstruction, the main factors can be described as follows. (1) In the beam experiment, the detector position itself is approximate, not very accurate, or some parts of the detector may deviate from the central position, and the deviated directions may be different, this will influence on the selection of coordinate system, and cause the error of position calculation of a certain part of the detector relative to the coordinate axis or coordinate origin. It requires to make the position calibration before the data processing, ensuring that the calculations are performed in the same coordinate system.
468
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
(2) Because of the material of the detector, when a particle hits the detector, many secondary particles may be produced, including some charged particles, we can not judge that the cluster is composed of what kinds of particles. Of course, it is possible that a cluster may be produced also by the interaction of a recoil particle with the detector. (3) In the data processing, the centroid position of the cluster is selected as the hit point, but the position that a particle actually hits the detector is not certainly the centroid position, there is a deviation existed. 250 200 150 100 50 0 -200
-100
0
100
200
300
400
500
Fig. 13 The reconstructed trajectory of an electron with a momentum of 100 GeV and an incident angle of 0◦
300 200 100 0 -100 -200 -300 -200
-100
0
100
200
300
400
500
Fig. 14 The reconstructed trajectory of a proton with a momentum of 400 GeV and an incident angle of 0◦
6.
SUMMARY
The track reconstruction has an important significance for the detection of dark matter particles, to recover the tracks of incident particles will help us to obtain the information
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
469
about the type, energy, intensity, and incident direction of incident particles, and to analyze further the generation of particles and the direction of the source. If these particles are produced by the decay or annihilation of dark matter, we can learn about the distribution of dark matter, as well as the correlation with the distribution of other kind matter. The most important points for the track reconstruction are: the deletion of pseudo-points, track search, track fitting, and the relevant corrections. Several key points should be taken into consideration in the process of track reconstruction, for example the detector performance (mainly the geometric configuration and accuracy (resolution)), energy loss, multiple scattering, etc., so that we can build a mathematic model to describe the particle track with a fairly good accuracy. In order to verify the accuracies of the model and algorithm, before analyzing the actual experimental data, the whole reconstruction process should be tested by using the simulated data, and the real values of estimated parameters in the simulated data should be known already, in favor of the test on the normalized residuals. The exact errors of the fitted parameters should be investigated carefully, which may be influenced by the pseudo-points or pseudo-tracks, as well as the method of their deletion, and the position error of the detector installation. This paper adopts mainly the energy-weighted least squares method and the straight line fitting method based on the Kalman filtering to reconstruct the incident particle track, by comparisons, it is found that the result (error) of track reconstruction is related with the energy of the incident particle. Generally, the higher the incident energy, the smaller the error of the result of track reconstruction. Further more, we find that in the high-energy region, the result obtained by the energy-weighted least squares linear fitting and the result obtained after the prediction, filtering and smoothing of the Kalman filter are basically equivalent, but in the low-energy region, the Kalman filtering has a relatively apparent superiority. The Kalman filter is a kind of powerful, flexible, and effective tool, it can be used for not only the track fitting, but also the calculations of optimum extrapolation and interpolation, the search and deletion of pseudo-points, the error adjustment and detector calibration using the actual data, as well as the merge of track segments. However, the Kalman filter has also its own defects, some times it is difficult to find a suitable initial value of Kalman filtering, if the initial angular resolution and momentum resolution of the detector are very poor, and the next detector is not very near, the first-step prediction may have the following problems: (1) It is very difficult to calculate precisely the amount of matter to be penetrated by the real particle track, because the predicted trajectory may be completely different from the real trajectory. (2) If the initial value deviates far from the real track, the linearity approximation of the track model will be invalid. (3) When the predicted track does not intersect with the surface of the next detector, the prediction procedure will be unable to perform.
470
LU Tong-suo et al. / Chinese Astronomy and Astrophysics 41 (2017) 455–470
A method to solve this kind of problems is to take the reference track derived by the pattern recognition method as the initial point of linear approximation. Namely, to calculate first the intersections of the reference track with all surfaces (the so-called reference state), then to make prediction, filtering, and smoothing by using the difference between the state and the reference state. Besides, the position measurement during the instrument installation is not perfectly accurate, the calibration of the instrument will have a certain influence on the result of track reconstruction, hence, in order to improve the effectiveness of the track reconstruction, it is also necessary to make calibration on the instrument. References 1
Chang J. ChJSS, 2014, 34, 550
2
Dong Y. F., Study of Data Acquistion System for the Silicon Detector on the Dark Matter Detection Satellite, Bejing: Institute of High-energy Physics of Chinese Academy of Sciences, 2015
3
Zhang Y. L., Study of Electromagnetic Calorimeter for Dark Matter Detection, Hefei: The University of Scence and Technology of China, 2011
4
Xie M. G., Guo J. H., Wu J., et al., AcASn, 2014, 55, 170
5
Xie M. G., Guo J. H., Wu J., et al., ChA&A, 2015, 39, 129
6
Guo J. H., Cai M. S., Hu Y. M., et al., AcASn, 2012, 53, 72
7
Zhang L., Guo J. H., Zhang Y. Q., AcASn, 2014, 55, 522
8
Zhang L., Guo J. H., Zhang Y. Q., ChA&A, 2015, 39, 380
9
Fruhwirth R., Regle M., High-energy Physics Data Analysis, Translatged by Zhu Y. S., Hefei: Press of the University of Science and Technology of China, 2011
10 11
Liu J., Mao Z. P., Li W. G., et al., High Energy Physics and Nuclear Physics, 1998, 22, 587 Wang X. L., Particle Detection Technique, Hefei: The University of Science and Technology of China, 2009
12
Mao Z. P,, Wang T. J., Wang S. Q., HEPNP, 1993, 17, 193
13
Grupen C., Shwartz B., Spieler H., Particle Detectors, Cambridge: Cambridge University Press, 2011
14
Kalman R. E., Transactions of the ASME-Journal of Basic Engineering, 1960, 35
15
Wang J. R., Instrument and Meter for Automation, 2013, 140
16
Gao F., Dong H. Q., System Engineering and Electronic Technology, 2006, 28, 775