-
We take the muon lifetime measurement experiment as a demonstration. Two plastic scintillators were used to build a coincidence detector, in order to distinguish cosmic muon events from ambient radiation events. Only particles with high enough energy are able to go through both layers of the plastic scintillators almost simultaneously. Since the energy of cosmic muons is
$ {\cal{O}}(1) \ \mathrm{GeV} $ , these particles will penetrate the scintillators with no difficulty and leave simultaneous coincident signals. However, the particles from ambient radiation ($ \alpha,\;\beta $ and$ \gamma $ ), which is one of the major detector noise sources, are of lower energy and will hardly be able to trigger both of the scintillators. The schematic diagram of the coincidence detector and the circuit diagram to readout SiPM signals for such a purpose are shown in Fig. 1 and Fig. 2. First, we will focus on the classification of the event types and their corresponding waveforms. Then, the statistical method will be applied to suppress the pile-up background in the experimental data analysis. -
To the best of our knowledge, there are mainly two kinds of muon events in the detector scintillators. The first one is the muon coincidence event, which is trivial and less interesting. The second one is the muon decay events. The cosmic muons of proper energy will stop inside the lower scintillator. Muons or antimuons are not stable and can decay into an electron or a positron and neutrino-antineutrino pairs. The decay process of
$ \mu^{\pm} $ is shown as follows:$$ \mu^{\pm} \rightarrow \mathrm{{e}}^{\pm} + \nu_e (\bar{\nu}_e) + \bar{\nu}_\mu (\nu_\mu). $$ (1) The secondary electron/positron can also ionize in the plastic scintillator and produce a second pulse signal. The time interval
$ \varDelta t $ between two pulse peaks of the arrival signal caused by a passing muon and the decay signal by Michel electron will shape an exponential distribution ideally as samples for the muon lifetime measurement.As for the backgrounds, we pay more attention to background sources from cosmic muon pile-up events, electronic noise and ambient radiation induced coincident events. The detector structure indicates that most of high-energy cosmic muons are likely to traverse the lower scintillator panel only. If the time interval between a coincidence event and a single-bottom-hit event is located outside the specified time window
$ T $ but fall within the sampling range in the readout system, these events are also saved and included in data analysis. However, we have to carefully deal with these pile up events to avoid contamination in the real cosmic muon decay signals. Another noise sources are$\alpha$ ,$\beta$ and$\gamma$ particles from the ambient radiation. The indoor ambient radioactivity level is set to 100$ \mathrm{{Bq/m}}^3 $ , higher than the average levels in China (~$ 40 \ \mathrm{{Bq/m}}^3 $ )[30-31]. Taking them into simulation, we immediately realize that these events are extremely rare in the current study and can be safely ignored. The responding mode of the coincidence events, the muon decay events, the muon pile-up events in shown in Fig. 3.Figure 3. (color online)The classification of the detector responses. The ''Event Type'' represents the actual type of an event. The normal muon coincidence events (Coin.), the muon decay events and the muon pile up events are all taken into account. The "Upper pulse" and ''Lower pulse" represent the pulse signals generated by the upper and lower scintillator, respectively. An OR logic gate is applied to two channels of each scintillator, in order to maximize the detection efficiency. The ''Event Classification" represents the reconstruction of the event type based on the waveform. The record time window is activated once the upper and lower SiPMs are almost simultaneously triggered and will keep recording for 20 μs. If there are 2 pulses within a single time window, the event will be classified as a muon-decay event.
-
For simplicity, we can start with a rough estimate on pile-up events in the cosmic muon lifetime measurement experiment. Assuming that the muon events are uniformly distributed in time, the peak distance of pile-up events, denoted as
$ D_{ \rm{pileup }}$ , will also almost uniformly distributes within the time window,$ T $ . Thus, we can estimate that the mean value of$ D_{ \rm{pileup}} $ satisfies$\bar{D}_{ \rm{pileup}}\approx \frac{1}{2}T > \tau_{\mu}$ . It means that misidentified pile-up events will lead to a larger value in the muon lifetime result.In fact, since the muon pile-up events can be treated as two adjacent muon coincidence events with small enough time interval, denoted as
$ \varDelta T $ . According to the previous assumption that the detector responses are uniformly distributed in time,$ \varDelta T $ obeys the exponential distribution and the probability distribution function (PDF) is as follows:$$ f_c^\infty (t) = \left\{ \begin{array}{*{20}{l}} \dfrac{1}{{{\tau _c}}}{{\rm{e}}^{ - t/{\tau _c}}}&(t \in [0,\infty ))\\ 0 &({\rm{elsewhere}}), \end{array} \right. $$ (2) where
$ \tau_c $ is arithmetic average of the$ \varDelta T $ samples and determined by data driven analysis.$$ {f_c}(t) = \left\{ \begin{array}{*{20}{l}} \dfrac{1}{{A{\tau _c}}}{{\rm{e}}^{ - t/{\tau _c}}}&(t \in [0,T])\\ 0 &({\rm{elsewhere}}), \end{array} \right. $$ (3) where
$ A = \int_{0}^{T} {f_c^{\infty}(t) \mathrm{d} t} = 1 - \mathrm{e}^{-T/\tau_c} $ , is the normalization factor. The muon lifetime by its definition obeys the exponential distribution. The two PDFs are shown in Fig. 4. According to our assumption, the pile-up events are almost uniformly distributed in the time window$ T $ , which means$ \tau_c $ , the mean value of$ \varDelta T $ satisfies$ \tau_c \gg T > \tau_\mu $ . As a result, the$ \varDelta T $ 's PDF is a slow-changing function in$ \left[0, T\right] $ and acts like a uniform distribution, as shown in the Fig. 4.Figure 4. (color online)Peak Distance's distributions of muon decay events and event pile up effects. The blue line and the red line represent the muon's lifetime distribution and the background events distribution in a given time window
$T = 13 \ \mu \text{s}$ , respectively. The gray area is the time window cut.Let random variables
$D,\ c$ and$ m $ represent the peak distance of double-pulse events, pile-up events, muon decay events, respectively. Without a loss of generality, we assume that the muon pile-up and decay events are independent. Allow$ \lambda $ to represent the proportion of muon decay events in all double-pulse events with$ \lambda \in \left(0, 1\right) $ . As a result, the peak distance of double-pulse events$ D $ can be expressed in the form of weighted sum of the peak distance of pile-up events and that of muon decay events, whose weight is$ 1-\lambda $ and$ \lambda $ , respectively, as shown in Eq. (4).$$ D = \lambda m + (1-\lambda) c \equiv M + C, $$ (4) where
$ M\equiv \lambda m, C\equiv (1-\lambda)c $ . Since the PDFs of$ m $ and$ c $ are clear, one can immediately write down the PDFs of$ M, \;C $ and$ D $ , as shown in Eqs. (5)~(7), respectively.$$ {f_M}(x) = \left\{ \begin{array}{*{20}{l}} \dfrac{1}{{{\tau^{} _1}}}{{\rm{e}}^{ - x/{\tau ^{}_1}}}&(x \in [0,\infty ))\\ 0 &({\rm{elsewhere}}), \end{array} \right. $$ (5) $$ {f_C}(y) = \left\{ \begin{array}{*{20}{l}} \dfrac{1}{{A{\tau _2^{}}}}{{\rm{e}}^{ - y/{\tau _2^{}}}}&(y \in [0,T])\\ 0 &({\rm{elsewhere}}), \end{array} \right.$$ (6) $$ \begin{split} {f_D}(z) =& \int\nolimits_{ - \infty }^{\,\infty} {{f_M}(x){f_D}(z - x){\rm{d}}x} \\ =& \left\{ \begin{array}{l} \dfrac{1}{{A({\tau _1^{}} - {\tau _2^{}})}}\left( {{{\rm{e}}^{ - \frac{z}{{{\tau _1}}}}} - {{\rm{e}}^{ - \frac{z}{{{\tau _2^{}}}}}}} \right)\;\;\left( {t \in \left[ {0,(1 - \lambda )T} \right)} \right),\\ \dfrac{1}{{A({\tau _1^{}} - {\tau _2})}}{{\rm{e}}^{ - z/{\tau _1^{}}}}\left\{ {1 - \exp \left[ { - \dfrac{{{\tau _1} - {\tau _2^{}}}}{{{\tau _1^{}}{\tau _2}}}(1 - \lambda )T} \right]} \right\}\\ \left( {t \in \left[ {(1 - \lambda )T,\infty } \right)} \right),\\ 0\;\;({\rm{elsewhere}}). \end{array} \right.\\[-13pt] \end{split} $$ (7) The mean value of D is shown in the Eq. (8).
$$ \bar{D} \approx E(D) = \int\nolimits_{0}^{\,\infty}{z f_D(z) \mathrm{d}z} = \tau_1^{} + \tau_2^{} - \frac{T(1-\lambda)}{A}\mathrm{e}^{-T / \tau_c} , $$ (8) where
$\tau_1^{} = \lambda \tau_\mu$ ,$ \tau_2 = (1-\lambda) \tau_c $ . Finally, we can replace the parameters$ A $ and$ \lambda $ with the observable quantities. If we allow$N, \;N_{\rm D},\; N_{\rm C}$ and$N_\mu$ to represent the number of coincidence events, double-pulse events, noise events, muon decay events, respectively, we will obtain the following equations:$$ N_{\rm C} = N \int\nolimits_{0}^{\,T} {f_c^{\infty}(t) \mathrm{d} t} = N A , $$ (9) $$ N_{\mu} = \lambda N_{\rm D}, $$ (10) $$ N_{\rm C} = (1-\lambda) N_{\rm D}. $$ (11) Using Eqs. (9)~(11) to replace the parameters A and
$ \lambda $ , we finally get to the result for$\tau_\mu$ :$$ \tau_{\mu} = \frac{\bar{D} + \dfrac{N}{N_{\rm D}}T\mathrm{e}^{-T/\tau_c} - \dfrac{N}{N_{\rm D}}\left( 1 - \mathrm{e}^{-T/\tau_c} \right) \tau_c}{1 - \dfrac{N}{N_{\rm D}}\left( 1 - \mathrm{e}^{-T/\tau_c} \right)}. $$ (12) Keep in mind that
$ \tau_c \gg T $ ,$ \frac{T}{\tau_c} \ll 1 $ . In order to interpret the physics in the above formulae, we can expand the term$ \mathrm{{e}}^{-T/\tau_c} $ in Taylor series. With the numerator to the first order and the denominator to the zeroth order, then one could reach:$$ \tau_{\mu} \approx \bar{D} - \frac{N}{N_{\rm D}}\frac{T^2}{\tau_c}. $$ (13) Eqs. (13) is simplified and relatively elegant. This approximation explicitly shows the outcome of background events represented by the term
$ \frac{N}{N_{\rm D}}\frac{T^2}{\tau_c} $ . This term also verifies our initial guess that the effect of pile-up events will lead to a larger value in the lifetime measurement result. It is then easy to understand how to suppress the pile up coincident backgrounds. -
As it is relatively rare to identify the decay muon events by coincident measurement of cosmic muons, it must be good to check the validity of this statistic correction method by means of large samples by Monte Carlo simulation. For this purpose, we make use of GEANT4 with modularized, ready-to-use physics lists with the class G4ModularPhysicsList[33-35]. The built-in class G4EmStandardPhysics and G4OpticalPhysics are used for Electro-magnetic processes and optical processes, respectively. The class G4StoppingPhysics is also used to simulate the
$ \mu^{-} $ capture, which is verified to be affecting the lifetime of$ \mu^{-} $ [36-37]. As for the background study, the pile up events and the ambient radiation events are both simulated. The ambient radiation events are extremely rare, probably due to the high threshold value and the low background radiation intensity in our experimental setup. -
The data used for generating cosmic muon spectrum is taken from CAPRICE 94 and CAPRICE 97 experiments[38], whose results have relatively small errors and are in good agreement with other experiments and spectrum shape equations[39-40]. The zenith angle distribution of cosmic muon is described by the shape of
$ \cos^{\alpha}\theta $ with$ \alpha \approx 2 $ , which has been verified by multiple measurements[41-42]. -
Geant4 can only track the semi-classic particles that have a well-defined position and momentum. In the simulation program, the generated light signal is a collection of single photoelectrons (SPEs). To reconstruct the actual electronic signal, a conversion algorithm is needed. We have built a model to simulate the waveform of SiPM, based on the single photoelectron response. First of all, a model converting the infinitesimal-width light signal is proposed based on the assumption that SiPM's response is the linear superposition of the SPEs. Second, we apply the model to convert the finite width light signal in the most general case.
The signal of SPE is fitted by a model with a convolution of logarithmic and exponential terms, as is shown in Eq. (14). For simplicity, we can first neglect the width of the light signal. Hence, it becomes the Dirac delta function,
$ N \delta (t) $ , where the parameter$ N $ represents the relative light intensity. Since the PN nodes inside a SiPM are in parallel with each other, it is reasonable to assume that if the PN nodes are ignited simultaneously, the total current is the linear superposition of each PN node's current. Therefore, the electronic signal can be writen down by a simple multiplication, as is shown in Eq. (15).$$ SPE(t) = \frac{A}{t\sigma \sqrt{2\pi }}\exp{\left[ -\frac{\ln^2 \left( \dfrac{t}{\tau} \right)}{2\sigma^2} \right]} + B, $$ (14) $$ y(t) = N\times SPE(t)\,, $$ (15) where
$ \tau $ and$ \sigma $ accounts for the central value and the width in the SPE waveform, respectively.$ A $ and$ B $ represents the amplitude and the baseline of each waveform, respectively. If the light signal has a finite width, which means that photons hit the SiPM at the different time, each hit can be treated as an individual infinitesimal-width light signal and trigger a SPE response. Hence, the reconstructed electronic signal is the linear superposition of all SPE responses. To reduce the computational resources, one can fill the time stamps of photon hits$ \left\{ t_i'\right\}_{i = 1}^{N} $ into a histogram and get the bin edges$ \left\{ t_i \right\}_{i = 1}^{n} $ and the corresponding bin counts$ \left\{ n_i \right\}_{i = 1}^{n} $ . Each bin can be treated as a infinitesimal-width light signal$ n_i \delta(t - t_i) $ . As a result, the light signal$ I(t) $ is the superposition of all bins and the reconstructed electronic signal$ Y(t) $ is the convolution of SPE waveform$ SPE(t) $ and the light signal$ I(t) $ , as shown in Eq. (16). Fig. 11 shows the typical reconstructed electronic signal.Figure 11. (color online)The simulated response of muon decay events. The reconstructed electronic signal is convoluted with the electronic noise.
$$ \begin{split} I(t) =& \sum\limits_{i = 1}^n {{n_i}} \delta (t - {t_i}),\\ Y(t) =& \int\nolimits_{ - \infty }^{\,\infty} {y(t)I(t){\rm{d}}t} = \sum\limits_{i = 1}^n {{n_i}} SPE(t - {t_i}). \end{split}$$ (16) -
In total, 33355 double-pulse events in simulation are collected and analyzed. The event pile up effect and the ambient radiation events are also included. After we apply the improved analysis method in the simulated data samples, the result with correction is
$ 2.20 \pm 0.03 $ µs at 95% C.L. while the simulation result without correction gives$ 2.29 \pm 0.03 $ µs. The results only include the statistical error. The histogram of muon decay spectrum in simulation is shown in Fig. 12. Therefore, the consistency of our result with regards to the reference value in simulation has verified the validity of the correction.
An Improved Treatment of Pile-up Events Demonstrated by a Cosmic Muon Lifetime Measurement Experiment
-
摘要: 提出了一种在数据分析中降低堆积事例误差的一种修正方法,并通过塑料闪烁体搭建的宇生μ子时间符合探测器开展寿命测量实验进行验证。寿命测量实验研究表明,主要本底事件来自电子学噪声和堆积事例。为了弥补在本地实验室中,μ子衰变事例在总体宇生μ子事例中较为稀有的短板,我们使用蒙特卡罗模拟程序产生大样本量,进一步验证该分析方法的有效性。修正后的最终实验结果为
$\tau_{\mu}^{\rm exp} = 2.19 \pm 0.07$ µs,而修正前的实验结果为$\tau_{\mu}' = 2.27 \pm 0.07$ µs。(95%置信度水平)。预计,统计学处理堆积事例的方法将适用于符合测量与堆积事例伴生的同类型实验。Abstract: A statistical correction method is proposed to suppress the pile-up background events in data analysis. This method is verified by a cosmic muon lifetime measurement experiment, achieved by a plastic scintillator detector using the timing coincidence method, where dominant background events originate from pile up muons and electronic noise. To complement the intrinsic shortcoming of relatively rare decay events from registered cosmic muon events in the local laboratory, Monte Carlo simulation is applied to generate large samples in order to cross check the new method. The measurement of the muon lifetime in our setup gives a result of$\tau_\mu^{\rm exp} = 2.19 \pm 0.07 \ \mu \text{s}$ at 95% confidence level, while the result before applying the correction is$\tau_{\mu}' = 2.27 \pm 0.07 ~\mu \text{s}$ (95% C.L.). The treatment of pile-up events by a statistical correction equation in this study might be adapted to improve data analysis in the general coincident background dominating experiments.-
Key words:
- cosmic muon /
- muon lifetime measurement /
- coincident backgrounds /
- pile-up events
-
Figure 3. (color online)The classification of the detector responses. The ''Event Type'' represents the actual type of an event. The normal muon coincidence events (Coin.), the muon decay events and the muon pile up events are all taken into account. The "Upper pulse" and ''Lower pulse" represent the pulse signals generated by the upper and lower scintillator, respectively. An OR logic gate is applied to two channels of each scintillator, in order to maximize the detection efficiency. The ''Event Classification" represents the reconstruction of the event type based on the waveform. The record time window is activated once the upper and lower SiPMs are almost simultaneously triggered and will keep recording for 20 μs. If there are 2 pulses within a single time window, the event will be classified as a muon-decay event.
Figure 4. (color online)Peak Distance's distributions of muon decay events and event pile up effects. The blue line and the red line represent the muon's lifetime distribution and the background events distribution in a given time window
$T = 13 \ \mu \text{s}$ , respectively. The gray area is the time window cut.Figure 6. (color online)The collected double-pulse waveforms from either end of the lower scintillator can be categorized into three sets. Set I waveforms only have a single main peak(MP) and a ringing peak(RP) (MP-RP waveforms), while both Set II and III have 2 main peaks (MP-MP waveforms). The STH represents the static threshold. A double-pulse wavefrom is recorded only if it has at least one peak that higher than STH, the static threshold.
Figure 7. (color online)The pulse height spectra of the first and the second peaks in the waveforms. The abscissa represent the first peaks'(mostly MP) height and the vertical one represents the second peak's (mostly RP) height. The width of the confidence bound is determined by the range at
$ 3\sigma$ C.L. The upper bound is used for the dynamic threshold, DTH. The excluded MP-MP events are represented by red circles.Figure 8. The flow chart to categorize Set I and III waveforms in Fig. 6. Set I waveforms have adjacent main peaks and ringing peaks, while Set III waveforms have two adjacent main peaks. Hence, Set I waveforms are also recorded and analyzed, though they are not decay signals. To distinguish Set I and III waveforms, the dynamic threshold DTH is introduced.
-
[1] BERTOLINI D, HARRIS P, LOW M, et al. JHEP, 2014, 10: 059. doi: 10.1007/JHEP10(2014)059 [2] AABOUD M, AAD G, ABBOTT B, et al. Eur Phys J C, 2017, 77(9): 580. doi: 10.1140/epjc/s10052-017-5081-5 [3] ARJONA MARTŃEZ J, CERRI O, PIERINI M, et al. Eur Phys J Plus, 2019, 134(7): 333. doi: 10.1140/epjp/i2019-12710-3 [4] SIRUNYAN A M, TUMASYAN A, ADAM W et al. JINST, 2020, 15(09): P09018. doi: 10.1088/1748-0221/15/09/P09018 [5] MAJUMDAR K. J Phys Conf Ser, 2015, 598(1): 012023. doi: 10.1088/1742-6596/598/1/012023 [6] ARMATOL A, ARMENGAUD E, ARMSTRONG W, et al. arXiv: 2011.11726, 2020. [7] SIRUNYAN A, TUMASYAN A, ADAM W, et al.Journal of Instrumentation, 2017, 12(10): P10003. [8] KOCH L, DENIG A, DREXLER P, et al. Journal of Instrumentation, 2017, 12(07): C07001. doi: 10.1088/1748-0221/12/07/c07001 [9] NAKHOSTIN M. Nuclear Engineering and Technology, 2020, 52 (2): 360. [10] LEE M, LEE D, KO E, et al. Nuclear Engineering and Technology, 2020, 52(5): 1029. [11] CASIMIRO E, ANJOS J C. Journal of Physics: Conference Series, 2008, 116: 012003. doi: 10.1088/1742-6596/116/1/012003 [12] CAO J, LUK K B. Nuclear Physics B, 2016, 908: 62. [13] AN F, BALANTEKIN A, BAND H, et al. Nucl Instr and Meth A, 2015, 773: 8. [14] SUZUKI A. The European Physical Journal C, 2014, 74(10). doi: 10.1140/epjc/s10052-014-3094-x [15] COLLABORATION R, AHN J K. arXiv: 1003.1391, 2010. [16] LUNDY R A. Phys Rev, 1962, 12: 1686. doi: 10.1103/PhysRev.125.1686 [17] BARDIN G, DUCLOS J, MAGNON A, et al. Phys Lett B, 1984, 137 (1): 135. [18] CHITWOOD D B, BANKS T I, BARNES M J, et al. Phys Rev Lett, 2007, 99: 032001. doi: 10.1103/PhysRevLett.99.032001 [19] ŚANCHEZ E. Physics Procedia, 2011, 17: 168. [20] SUN Lazhen, WU Yusheng, LI Cheng. Physics Experiment Ation, 2010, 30(2): 1. (in Chinese) doi: 10.3969/j.issn.1005-4642.2010.02.001 [21] HALL R E, LIND D A, RISTINEN R A. American Journal of Physics, 1970, 38(10): 1196. doi: 10.1119/1.1976002 [22] OWENS A, MACGREGOR A E. American Journal of Physics, 1978, 46(8): 859. doi: 10.1119/1.11407 [23] LEWIS R J. American Journal of Physics, 1982, 50(10): 894. doi: 10.1119/1.13013 [24] COAN T, LIU T, YE J. American Journal of Physics, 2006, 74(2): 161. doi: 10.1119/1.2135319 [25] WEBBER D M, TISHCHENKO V, PENG Q, et al. Phys Rev Lett, 2011, 106: 041803. doi: 10.1103/PhysRevLett.106.041803 [26] WANG Rong, HAN Chengdong, ZHANG Yapeng, et al. Nuclear Physics Review, 2016, 33(03): 315. (in Chinese) [27] TIAN Yi, HU Luguo, SUN Baohua. College Physics, 2018, 37(10): 36. (in Chinese) [28] LIN Yanchang, CHEN Shaomin, GAO Yuanning, et al. Experimental Technology and Management, 2008, 25(09): 41. (in Chinese) [29] WU Yusheng, LV Zhiyan, LI Shu, et al. Journal of University of Science and Technology of China, 2010, 040(006): 608. (in Chinese) [30] ORGANIZATION W H. Who Handbook on Indoor Radon: a Public Health Perspective[M]. Switzerland: WHO Press, 2009: 89. [31] PAN Z. Radiation Level in China[M]. Beijing: China Atomic Energy Publishing and Media Co., Ltd, 2011. [32] GROUP P D, ZYLA P A, BARNETT R M, et al. Progress of Theoretical and Experimental Physics, 2020, 2020(8). doi: 10.1093/ptep/ptaa104 [33] AGOSTINELLI S, ALLISON J, AMAKO K, et al. Nucl Instr and Meth A, 2003, 506(3): 250. [34] ALLISON J, AMAKO K, APOSTOLAKIS J, et al. IEEE Transactions on Nuclear Science, 2006, 53(1): 270. doi: 10.1109/TNS.2006.869826 [35] ALLISON J, AMAKO K, APOSTOLAKIS J, et al. Nucl Instr and Meth A, 2016, 835: 186. [36] BARDIN G, DUCLOS J, MAGNON A, et al. Nuclear Physics A, 1981, 352(3): 365. [37] WARD T, BARKER M, BREEDEN J, et al. American Journal of Physics, 1985, 53(6): 542. doi: 10.1119/1.14235 [38] KREMER J, BOEZIO M, AMBRIOLA M L, et al. Phys Rev Lett, 1999, 83: 4241. doi: 10.1103/PhysRevLett.83.4241 [39] GAISSER T K, ENGEL R, RESCONI E. Cosmicrays and Particle Physics[M]. Cambridg: Cambridge University Press, 2016. [40] GUAN M, CHU M C, CAO J, et al. arXiv: 1509.06176, 2015. [41] BAHMANABADI M. Nucl Instr and Meth A, 2019, 916: 1. [42] PETHURAJ S, DATAR V, MAJUMDER G, et al. Journal of Cosmology and Astroparticle Physics, 2017, 2017(09): 021. doi: 10.1088/1475-7516/2017/09/021