aa r X i v : . [ a s t r o - ph . H E ] A p r Chinese Physics C Vol. XX, No. X (2019) XXXXXX
The expectation of cosmic ray proton and helium energyspectrum below 4 PeV measured by LHAASO
L.Q. Yin , S.S. Zhang , Z. Cao , , B.Y. Bi , C. Wang , J.L. Liu , L.L. Ma ,M.J. Yang , Tiina Suomij¨arvi , Y. Zhang , Z.Y. You , , Z.Z Zong for the LHAASO Collaboration Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,100049 Beijing, China. University of Chinese Academy of Science, Beijing, China. Institut f¨ur Astronomie und Astrophysik, Eberhard Karls Universit¨at, T¨ubingen 72076, Germany. Kunming University, Kunming, Yunnan, China. Institut de Physique Nucl´eaire d’Orsay, IN2P3-CNRS, Universit´e Paris-Sud, Universit´e, Paris-Saclay, 91406Orsay Cedex, France. China Nuclear Power Engineering Co., Ltd., 100084, Beijing, China
Abstract:
Large High Altitude Air Shower Observatory(LHAASO) is a composite cosmic rayobservatory consisting of three detector arrays: kilometer square array (KM2A) which includesthe electromagnetic detector array and muon detector array, water Cherenkov detector array(WCDA) and wide field of view Cherenkov telescope array (WFCTA). One of the main scientificobjectives of LHAASO is to precisely measure the cosmic rays energy spectrum of individualcomponents from 10 eV to 10 eV. The hybrid observation will be employed by LHAASOexperiment, in which the lateral and longitudinal distributions of the extensive air showercan be observed simultaneously. Thus many kinds of parameters can be used for primarynuclei identification. In this paper, high purity cosmic ray simulation samples of light nucleicomponent are obtained through Multi-Variable Analysis. The apertures of 1/4 LHAASO arrayfor pure proton and mixed proton and helium (H&He) samples are 900 m Sr and 1800 m Sr respectively. A prospect of proton and H&He spectra from 100 TeV to 4 PeV is discussed. Key words:
LHAASO; hybrid measurement; energy spectrum; composition; TMVA.
An unsolved problem in cosmic ray observation is the “knee” structure in the energy spectrum,namely a significant bending of the spectrum from the power-law index of approximately -2.7to -3.1 around a few PeV. The causes of the knee structure are closely related to the origin,acceleration and propagation mechanism of cosmic rays, but there is no consistent measurementresult so far. Several experiments have measured the all-particle spectra around the knee region,such as KASCADE-Grand [1], Tibet-AS γ [2], EAS-TOP [3], etc. These experiments exhibited asimilar knee-like structure at energies of about 4 PeV, within a factor of two in the flux values [4].In addition, the energy spectra of individual components are much different. For example,the unfolded proton spectrum measured by the KASCADE experiment shows the steepening at2 PeV or 3PeV based on QGSJET01 model and SIBYLL2.1 model respectively [5]. However,the hybrid experiment of ARGO-YBJ and two prototypes of WFCTA shows that the energyspectrum for H&He has a knee-like structure at 700 TeV [6]. The main reasons for this situationare the lack of absolute energy scale and the way to identify the type of the primary particles ofcosmic rays. LHAASO [7, 8], located in Daocheng Haizishan, 4410 m a.s.l., Sichuan Province,China, will throw light on this traditional problem.
1) E-mail: [email protected]
XXXXXX-1hinese Physics C Vol. XX, No. X (2019) XXXXXX
The LHAASO site is at an ideal altitude for knee physics research because the atmosphericdepth is close to the maximum of the development of cosmic ray extensive air shower (EAS) [9],so that the fluctuation of EAS and the dependence of the interaction model will be the minimum.Similar to the ARGO-YBJ experiment, WCDA can obtain absolute-energy-scale through moonshadow observation [10]. Because the detection threshold of WCDA can be as low as severalhundred GeV, it can be directly compared with the measurement results of space experimentsto study the error of energy scale.Furthermore, LHAASO contains four types of detectors, which can detect electromagneticparticles, muons, Cherenkov and fluorescent photons in the EAS. Multi-Variable measurementsof the EAS can effectively distinguish the components of cosmic rays. Thus, the energy spectraof the individual components can be precisely measured by LHAASO through Multi-VariableAnalysis (MVA).LHAASO is planed to obtain the consecutive measurement of energy spectra by four stages.First of all, WCDA delivers the absolute-energy-scale to WFCTA by hybrid observation of cosmicrays from 10 TeV to 100 TeV. In this energy band, the energy spectra measured by LHAASOcan overlap with the spectra of space experiments. The second stage is mainly for the kneeobservation of light nuclei components of cosmic rays in the energies from 100 TeV to 10 PeV bythe hybrid detection of WFCTA, WCDA, and KM2A. In the third stage, the layout of WFCTAwill be changed to measure the knee of the heavy nuclei from 10 PeV to 100 PeV [11]. WCDAwill not be included in the hybrid detection due to saturation. The last stage is for the secondknee study from 100 PeV to 1 EeV, in which WFCTA will be operated in fluorescence observationmode [12].In the first three observation stages, the elevation angles of the WFCTA are 90 o , 60 o , 45 o respectively. Their corresponding atmospheric depths are respectively around the shower max-imum positions of cosmic rays in their own observation stages. This paper is devoted to theprospect of energy spectra observation of pure proton and H&He in the second stage throughsimulation. LHAASO was formally approved on December 31, 2015 and now is under construction. Thedetectors layout of LHAASO experiment is shown in Fig. 1. LHAASO consists of four types ofdetectors: electromagnetic detectors (red dots), muon detectors (blue dots), water Cherenkovdetector (blue squares) and wide field of view (FOV) Cherenkov telescopes(black rectangles).There are many vacancies in the muon detectors (blue dots) distribution, which is due to theterrain making it impossible to install detectors.WCDA [13], located at the center of LHAASO, is a water Cherenkov detector array with atotal area of 78,000 m . It detects the Cherenkov light produced in water by cascade processes ofsecondary particles in EAS. WCDA consists of three water ponds: two have area of 150 m × m and one has area of 300 m × m . All of the ponds have the water depth of 4.5 m. Each pondis divided into small cells with an area of 5 m × m , having two Photomultiplier Tubes (PMT)anchored on the bottom at the center of the cell. In one of the 150 m × m ponds, the pondin the lower left quarter in Fig. 1, a 1-inch PMT is placed right by an 8-inch PMT to enlargethe dynamic range to 200,000 photoelectrons. It allows the measurement of cosmic rays withenergy up to 10 PeV [14]. This pond is called WCDA++ to achieve the hybrid detection withthe WFCTA.WFCTA is located at 85 meters away from the center of WCDA++. It detects the Cherenkovand fluorescence photons in the air shower. The angle of elevation of the telescope’s main axis isset to 60 o . Each telescope includes a spherical mirror to collect Cherenkov and fluorescence pho- XXXXXX-2hinese Physics C Vol. XX, No. X (2019) XXXXXXFig.
1: The detectors layout of LHAASO experiment.tons and reflect them onto the camera [15] with the effective area of 5 m . The camera is locatedat the focal point of the spherical mirror and its function is to convert Cherenkov/fluorescencephotons into electrical signals. The optical sensor of the camera is SiPM [16], with 1024 piecesin total. The pixel of each SiPM is 0 . ◦ × . ◦ and the camera watches a FOV of 16 ◦ × ◦ .All the parts of the telescope are placed in a container for the convenience of movement andarrangement. Two telescope prototypes have been operated successfully at Yangbajing cosmicray observatory in Tibet [17].KM2A contains two sub-arrays: electromagnetic detector array (ED) and muon detector array(MD). ED array [18] covers an area of 1.3 km consisted of 5195 plastic scintillator detectors witha spacing of 15 meters. One ED detector has an effective area of 1 m . As its name describes,the ED detects the electromagnetic particles in the EAS. MD array is an underground waterCherenkov detector array [19]. There are 1171 muon detectors also covering 1 km area with aspacing of 30 meters. The area of one muon detector is 36 m , and it detects the information ofmuon in the EAS.The 1/4 LHAASO array includes the WCDA++, six Cherenkov telescopes, 1272 EDs and300 MDs. The layout of the 1/4 array is shown in Fig. 1. This 1/4 array will run for a few yearsto achieve the physical goal of light nuclei components spectra measurement.Since full-coverage array can provide more accurate geometric information of air showers, thesecondary particles are mainly measured by WCDA++ in the second hybrid detection stage. Itmeans that the shower core position in this stage is outside of the ED array. Therefore, ED isnot included in this simulation. The cascade processes of primary cosmic rays in the atmosphere are simulated by the COR-SIKA [20] program with the version of 6990. The EGS4 model is chosen for the electromagneticinteraction. The QGSJET02 and GHEISHA models are chosen for the high and low energyhadronic processes respectively. Both the information of Cherenkov photons and secondary par-ticle at the level of observatory are recorded to simulate the hybrid observation. Five components,proton, helium, CNO, MgAlSi, and iron are generated with energies from 10 TeV to 10 PeV ac-cording to a power law spectrum with an index of -2.7. The directions of the showers are from 24 ◦ XXXXXX-3hinese Physics C Vol. XX, No. X (2019) XXXXXX to 38 ◦ in zenith and from 77 ◦ to 103 ◦ in azimuth. The shower core position is evenly distributedin an area of 260 m ×
260 m. The primary energy spectrum is normalized to the exponent of -2.7from 10 TeV to 10 PeV, as shown in Fig. 2. Different components are normalized to the samewith proton.In addition, five components are also normalized according to H¨orandel model. The relation-ship between the original energy spectra and the two normalized spectra is that the number oforiginal proton events, the number of proton events with exponent of -2.7 and the number ofproton events with exponent from H¨orandel model are the same at 100 TeV. The events num-ber of other energy bins and other components are normalized in proportion. The normalizedstatistics of each component are equivalent to the observation data of one year with 10% dutycycle according to H¨orandel model [21].The simulation of the response of three LHAASO detector arrays for EAS is performedseparately. In the simulation, the actual coordinates of the detector arrays are transferred tothe CORSIKA coordinate system. The program of photons ray-tracing is used for WFCTAsimulation. The trigger pattern of the telescope is set to 7, namely seven neighboring pixelsshould be triggered. The fast simulation program [22] is used for the simulation of WCDA++and MD. Only photoelectrons are sampled according to the energy, direction and particle id ofthe secondary particles in EAS. The time response is not included. Then the simulated resultsof different detectors are merged before reconstruction. (E/GeV) log N u m be r o f E v en t s
10 protonheliumCNOMgAlSiiron
Normalized Spectrum
Fig.
2: The energy distribution of all simulation events after normalization.
According to the actual data acquisition, as long as the Cherenkov telescope is triggered, theevents measured by three detector arrays are reconstructed.Firstly, the shower core position is given by WCDA++ through NKG fitting [23]. Due to thelack of time response in WCDA++ fast simulation program, there is no information about theshower fronts. Therefore the arrival direction of the shower is obtained by a Gauss sampling, ofwhich angular resolution is 0.3 ◦ .Secondly, the perpendicular distance between the telescope and the shower axis ( R p ) is cal-culated. Next is Cherenkov image cleaning. The first step is to remove the fired pixels whosephotoelectrons measured is less than 30 in the image. The second step is to find the fired pixelwith the largest number of photoelectrons, and use it as the center to traverse the whole image XXXXXX-4hinese Physics C Vol. XX, No. X (2019) XXXXXX outward. Then remove all the islanded fired pixels.After image cleaning, the Hillas parameters [24], are given, as well as the total photoelectrons( N pe ) in the image. N pe is a good energy estimator. Before energy reconstruction, N pe shouldbe normalized to R p = 0 and α = 0, namely N pe = lg ( N pe ) + a × ( R p / m ) + b × tan ( α ) (1)where α is the space angle between the shower direction and the optical axis of the telescope,and the parameters a and b depend on the primary component of cosmic rays.Finally, the lateral distribution of muons is fitted by an empirical formula Eq.( 2) throughmaximum likelihood fitting [9]. The normalized muon size is given. ρ ( r, N µ ) = k G N µ ( rr G ) − k (1 + rr G ) − k [ m − ] (2)where k ∼ = 1 . k ∼ = 1 . r G ∼ = 220 m . These three parameters are reconfirmed according to thelayout of MD array. There are three basic principles for events selection: firstly, the shower core position falls inWCDA++; secondly, the Cherenkov image is complete; thirdly, the muon lateral distribution iswell fitted.In detail, the reconstructed shower core position located in 130 m × m arround WCDA++is selected. For a Cherenkov image, the number of fired pixels are more than 10; and the angulardistance between the weighted center of the image and the center of the camera plane is lessthan 6 ◦ . For muon lateral distribution fitting, the fitted normalized muon size ( k G N µ ) should bemore than 10 − .The number of events in each energy bin after selection is shown in Fig. 3 (left). The dottedline represents the number of proton events involved in the simulation, which is consistent withthe black line in Fig. 2. The solid lines represent the number of events after selection, and differentcolors represent different components. It points out that the hybrid observation becomes nearlyfully efficient above 100 TeV for all components. (E/TeV) log N u m be r o f E v en t s -1 protonheliumCNOMgAlSiiron proton before selection (E/TeV) log D e t e c t E ff ( % ) H & Heproton
Fig.
3: Left: the energy distribution after event selection. Right: the detection efficiency ofproton and H&He.Detection efficiency is defined as the ratio of the number of selected events to the number ofsimulation events in each energy bin. That is the black solid line divided by the blue dotted line
XXXXXX-5hinese Physics C Vol. XX, No. X (2019) XXXXXX in Fig. 3 (left). In this simulation, the detection efficiency is about 15% above 100 TeV for allcomponents, as shown in Fig. 3(right). The red dots show the detection efficiency of H&He andthe black dots show the the detection efficiency of proton. Besides, it is obvious that the errorbars of the last two points in Fig. 3 (right) are large due to the small statistics. Therefore theresults of lg ( E/T eV ) > . In this section, the parameters of the three types of detector arrays are analyzed according tothe development characteristics of the EAS. The component sensitive parameters of LHAASOhybrid observation are introduced. Then, the Toolkit for Multivariate Analysis (TMVA) [25, 26]is used and the selection results are presented. X max , is a traditional mass sensitive parameter. Here, ∆ θ , the exactly parameter usedto reconstruct the X max , is applied instead of the reconstructed X max . ∆ θ is the angular distancebetween the shower arriving direction and the gravity center of the Cherenkov image. But itcan not be used directly to classify primary particles because of the R p and shower energy ( N pe )dependence. After normalization, the structure of the mass sensitive parameter P X is as follows: P X = ∆ θ − . × R p − . × N pe (3)Moreover, the proton-induced showers exhibit an elongated elliptic shape [24]. And the ratioof length and width of the cherenkov image is also a traditional and effective parameter: P C = L/W − . × R p + 0 . × N pe (4)The distributions of P X and P C for proton and iron showers are shown in Fig. 4. X P -4 -3 -2 -1 0 1 N u m be r o f E v en t s protoniron Mean: -1.25RMS: 0.41 proton
Mean: -1.71RMS: 0.28 iron C P -2 -1 0 1 2 3 4 5 6 N u m be r o f E v en t s protoniron Mean: 2.81RMS: 0.48 proton
Mean: 1.84RMS: 0.37 iron
Fig.
4: The distributions of mass sensitive parameters P X (left) and P C (right) for proton (redline) and iron (blue) initialed showers.4.1.2 Parameters from MDThe muon size N µ of a shower heavily depends on the atomic number ( A ) of the primaryparticle: N Aµ /N pµ ≈ A (1 − η ) , where η is approximately 0.9 [9] and p indicates the proton shower. XXXXXX-6hinese Physics C Vol. XX, No. X (2019) XXXXXX
Therefore, the muon size by fitting ( N Fµ ), the total number of detected muons ( N Mµ ) and thenumber of fired MDs (NMD) are significant variables to identify the primary particles: P µ = lg ( N Fµ ) − . × N pe (5) P µ = lg ( N Mµ ) + 0 . × Rlp − . × N pe (6) P µ = lg ( N M D + 0 . × Rlp ) − . × N pe (7) Rlp is the perpendicular distance between the center of the MD array and the shower axis. Thedistributions of three parameters P µ , P µ and P µ for proton and iron showers are shown inFig. 5.4.1.3 Parameters from WCDAIt is also known that the iron initialed shower is more extensive in the same observation level.It means that the lateral extension of secondary particles in iron-induced air shower is larger.WCDA++, the fully covered detector array, can well detect the lateral distribution of secondaryparticles. The formulas for the average lateral distribution < ER > and fluctuation RM S areexpressed as Eq. 2 to Eq. 5. < ER > = P R i × P e i P P e i (8)where R i is the distance between the shower core position and the ith fired cell in WCDA++. P e i is photoelectronic measured by the ith fired cell. RM S x = s P x i × P e i P P e i − x weight (9) RM S y = s P y i × P e i P P e i − y weight (10) RM S = s P R i × P e i P P e i − < ER > (11)where the x weight and y weight are the centroid of distribution of secondary particles in WCDA++: x weight = P x i × P e i P P e i (12) y weight = P y i × P e i P P e i (13)These parameters are not only related to the primary component but also to the shower directionand the shower energy, hence it should be corrected before MVA: P F = lg < ER > − . × θ rec + 0 . × N pe (14) P F = RM S x − . × θ rec + 2 . × N pe (15) P F = RM S − . × θ rec + 1 . × N pe (16)The distributions of three parameters P F , P F and P F for proton and iron showers areshown in Fig. 5. XXXXXX-7hinese Physics C Vol. XX, No. X (2019) XXXXXX µ P -9 -8.5 -8 -7.5 -7 -6.5 -6 N u m be r o f E v en t s protoniron Mean: -7.41RMS: 0.23 proton
Mean: -6.90RMS: 0.19 iron µ P -4 -3.5 -3 -2.5 -2 -1.5 -1 N u m be r o f E v en t s protoniron Mean: -2.63RMS: 0.19 proton
Mean: -2.17RMS: 0.15 iron µ P -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 N u m be r o f E v en t s protoniron Mean: -1.16RMS: 0.11 proton
Mean: -0.87RMS: 0.10 iron F2 P N u m be r o f E v en t s protoniron Mean: 1.60RMS: 0.11 proton
Mean: 1.82RMS: 0.05 iron F3 P N u m be r o f E v en t s protoniron Mean: 27.27RMS: 2.88 proton
Mean: 34.93RMS: 2.39 iron F4 P
15 20 25 30 35 N u m be r o f E v en t s protoniron Mean: 24.99RMS: 1.89 proton
Mean: 27.95RMS: 0.97 iron
Fig.
5: The distributions of mass sensitive parameters P µ (up-left), P µ (up-right), P µ (middle-left), P F (middle-right), P F (bottom-left) and P F (bottom-right) for proton (red line) and iron(blue) initialed showers. XXXXXX-8hinese Physics C Vol. XX, No. X (2019) XXXXXX
In addition, nearby the shower core axis, the particle density of iron-shower is much less thanthat of proton-shower [6]. Hence the photoelectrons in the brightest cell ( N max ) measured byWCDA++ is sensitive to components: P F = lg ( N max ) − . × N pe (17)The total photoelectrons ( N pew ) measured by WCDA++ is also a mass sensitive variable, even itis weaker than N max : P E = lg ( N pew ) − . × N pe (18)The distributions of two parameters P F and P E for proton and iron showers are shown in Fig. 6. F P -5 -4 -3 -2 -1 0 1 N u m be r o f E v en t s protoniron Mean: -1.58RMS: 0.49 proton
Mean: -2.27RMS: 0.33 iron E P -2 -1.5 -1 -0.5 0 0.5 1 1.5 N u m be r o f E v en t s protoniron Mean: 0.24RMS: 0.29 proton
Mean: -0.02RMS: 0.18 iron
Fig.
6: The distributions of mass sensitive parameters P F (left) and P E (right) for proton (redline) and iron (blue) initialed showers.Other parameters have also been studied, such as the photoelectrons measured by the bright-est pixel in air Cherenkov image ( P C ), the total photoelectrons located 45 meters away fromthe shower core in the WCDA++ ( P F ), the muon density detected from 80 m to 100 m awayfrom the shower core position ( P µ ), and etc. Nevertheless, the particle classification of thesevariables is weak. After tuning of parameters, these ten parameters are candidates used as inputparameters for the BDTG classifiers.4.1.4 Correlation of ParametersThe parameters described above are not independent. There is a correlation between someparameters. Parameters provided by the same detector array are highly correlated, such as, P F and P F , P µ and P µ . The correlations of parameters from WCDA and MD for five componentsof cosmic rays are shown in Fig. 7. The parameters of P F and P F are negative correlation (leftplot) and P µ and P µ are position correlation (right plot). According to the development char-acteristics of the EAS, it is easy to understand that the more extensive of the lateral distributionof the secondary particles, the less particles nearby the shower core axis; and the more muondetectors fired, the more muon contents in the shower.The correlation of parameters provided by the different detector arrays is weak, such as, P F and P X , P µ and P F . The correlations of parameters from different detector arrays forfive components are shown in Fig. 8. The distributions are approximate to circle. However,the correlation between parameters of different components is slightly different; parameters forproton tend to be more correlated, as shown by the black dots in Fig. 8.The correlation matrix of these ten variables for proton initialed showers is shown in Fig. 9.The linear correlation coefficient 100% means a linear positive correlation and -100% means linear XXXXXX-9hinese Physics C Vol. XX, No. X (2019) XXXXXX F2 P F P -3.5-3-2.5-2-1.5-1-0.500.51 protonheliumCNOMgAlSiiron µ P -9 -8.5 -8 -7.5 -7 -6.5 -6 µ P -1.6-1.4-1.2-1-0.8-0.6-0.4-0.2 protonheliumCNOMgAlSiiron Fig.
7: The correlations of mass sensitive parameters from WCDA (left) and MD (right). X P -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 F P -4-3-2-101 protonheliumCNOMgAlSiiron F4 P
10 15 20 25 30 35 µ P -8.5-8-7.5-7-6.5-6 protonheliumCNOMgAlSiiron Fig.
8: The correlations of mass sensitive parameters from different detector arrays.negative correlation. 0 means no correlations among two parameters. It should be noticed thatthe parameter P C is basically independent of other parameters, as shown in the third column orin the eighth row in Fig. 9. This is because P C reflects the overall development characteristic,not just the lateral or longitude distribution of the air shower.After removing the highly correlated (97% or 98% correlations) parameters, six parameters, P C , P X , P F , P F , P µ and P µ , are used for TMVA training and analysis. Drawing on the experience of ARGO-YBJ & WFCT hybrid detection, two parameters areused for particle identification [27] in the event-by-event cut. However, too many variables areno longer suitable for particle identification through event-by-event cut because it is complicatedto get the optimal cut value of each parameter. Therefore, MVA is taken.As an important branch of statistics, multivariate analysis has been applied to most of thedisciplines. TMVA is specially developed for high energy physics based on the ROOT-integratedenvironment. It is powerful for signal and background classification. In accelerator physics, itcan effectively screen out the b-tagging signals from a large number of background particles inthe jet [28]. And likewise, it enables to identify the components of cosmic rays [29].TMVA works based on machine learning. It integrates multiple advanced algorithm classi-fiers, such as Boosted Decision Trees (BDT), Artificial Neural Networks (ANN), Support VectorMachine (SVM), and etc. Users need to input variable samples and select the algorithm. Thenthe machine training and testing are carried out. The result is to output one variable for usersto select signals.Here Boosted Decision Trees with Gradient boosting(BDTG) is chosen, which is the most
XXXXXX-10hinese Physics C Vol. XX, No. X (2019) XXXXXX -100-80-60-40-20020406080100 F P E P C P X P µ P µ P µ P F P F P F P F P E P C P X P µ P µ P µ P F2 P F3 P F4 P Correlation Matrix (signal)
100 97 10 72 59 58 50 -93 -92 -92 97 100 7 75 65 64 56 -91 -90 -91 10 7 100 31 -5 -3 -3 -17 -18 -15 72 75 31 100 50 52 39 -70 -68 -70 59 65 -5 50 100 97 91 -47 -44 -46 58 64 -3 52 97 100 93 -47 -44 -45 50 56 -3 39 91 93 100 -40 -37 -37-93 -91 -17 -70 -47 -47 -40 100 97 98-92 -90 -18 -68 -44 -44 -37 97 100 97-92 -91 -15 -70 -46 -45 -37 98 97 100
Linear correlation coefficients in %
Fig.
9: Linear correlation matrix for the input variables of proton events.widely used so far. The classification of H and H & He from other mass components is carriedout independently.The selected hybrid events are divided into two equal parts. One is used for TMVA trainingand testing; the other is used as data. The statistics of signal (proton or H&He) and background(others) are shown in the Table 1.Table 1: The number of events for BDTG classifier training and testing.NO. of events Proton H&HeSignal 40215 61098Background 91330 70447To avoid overtraining, several parameters are adjusted to achieve the best performance of theclassifier. Parameters of BDTG classifiers for the separation of H/other nuclei are as follows:- Number of trees in the forest is 380;- Minimum training events required in a leaf node are 30;- Max depth of the decision tree is 2.Parameters of BDTG classifiers for the separation of H&He/other nuclei are as follow:- Number of trees in the forest is 300;- Minimum training events required in a leaf node are 50;- Max depth of the decision tree is 2.The other parameters are set as default values.The training results are shown in Fig. 10. H & He can be well separated from other compo-nents; However, separation of proton from other nuclei is barely satisfactory. The backgroundrejection versus signal efficiency (“ROC curve”) is obtained, by cutting on the BDTG classifier XXXXXX-11hinese Physics C Vol. XX, No. X (2019) XXXXXX outputs for the events of the proton and H&He samples, as shown in Fig. 11. It is obviouslythat under the same signal efficiency, the background rejection of other heavy nuclei (H&He vs.other nuclei) is higher.
BDTG response -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 d x / ( / N ) d N Signal (test sample)Background (test sample) Signal (training sample)Background (training sample)
Kolmogorov-Smirnov test: signal (background) probability = 0.049 ( 0.07) U / O -f l o w ( S , B ): ( . , . ) % / ( . , . ) % TMVA overtraining check for classifier: BDTG BDTG response -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 d x / ( / N ) d N Signal (test sample)Background (test sample) Signal (training sample)Background (training sample)
Kolmogorov-Smirnov test: signal (background) probability = 0.934 (0.225) U / O -f l o w ( S , B ): ( . , . ) % / ( . , . ) % TMVA overtraining check for classifier: BDTG
Fig.
10: Training results of BDTG classifier for proton (left) and H&He (right).
Signal efficiency B a ck g r ound r e j e c t i on H & HeProton
Fig.
11: ROC curves BDTG classifier for proton (red line) and H&He (black line).
The training results are applied to the other half of the data sample. Generally, the signalevents can be selected according to the best cut points provided by TMVA. However, becauseour final goal is to obtain the selection efficiency under different energy bins, the distribution ofthe output of BDTG classifier versus reconstructed energy is studied firstly. As shown in Fig. 12,the separation of the output BDTG between H&He and heavy nuclei become larger as primaryenergy increases. The mean value of the signal (red dots) in each energy bin is slightly awayfrom the background (black dots) and the RMS of the signal is gradually decreased. The causesfor the separation becoming larger is that the performance of the input parameters gets betterin higher energies.To get a coincident selection efficiency, the cut values are scaled following the Eq. 19 andEq. 20:
CutV alue ( H ) = 0 . × lg ( Energy/T eV ) + 0 .
55 (19)
CutV alue ( H & He ) = 0 . × lg ( Energy/T eV ) − .
26 (20)
XXXXXX-12hinese Physics C Vol. XX, No. X (2019) XXXXXX /TeV)
Rec (Energy log B D T G -2-1.5-1-0.500.511.522.5 H & HeHeavy Nuclide
Fig.
12: The distribution of output parameter BDTG as a function of reconstructed energy. Theerror bars in the figure shows the RMS in each bin.After the cut described above, the aperture and contamination of hybrid observation arecalculated. The aperture of pure proton is about 900 m Sr and the aperture of H&He is about1800 m Sr , as shown in Fig. 13 (left). The contamination of proton sample is about 10% andthe contamination of H&He sample is less than 3% according to the H¨orandel model, as shownin Fig. 13 (right). (E/TeV) log s r) • A pe r t u r e ( m ProtonH & HeTMVA-BDTG (E/TeV) log C on t a m i na t i on ( % ) ProtonH & HeTMVA-BDTG
Fig.
13: The final selected aperture (left) and contamination (right) for H and H&He.Considering the uncertainty due to the hadronic interaction models, there are two batches ofdata, QGSJET-FLUKA and EPOS-FLUKA, were used for error analysis. It is found that theaperture of H&He is consistent within ± After the primary particle identification, the accuracy of event reconstruction for light com-ponents is obtained. For proton and H&He, the shower core resolution is less than 3 m. Theenergy reconstruction is obtained by WFCTA, as described in Section 3. There is a linear cor-relation between shower energy E and the corrected number of Cherenkov photoelectrons, N pe .For proton, a = 0 . b = 0 . E rec = 10 . × N pe − . . Thevariables a and b for H&He energy reconstruction are approximately equal to that of proton.The reconstructed bias and resolution for proton and H&He are shown in Fig. 14. The Bias isabout 4% and the resolution is less than 20%. XXXXXX-13hinese Physics C Vol. XX, No. X (2019) XXXXXX (E/GeV) log B i a s and R e s o l u t i on ( % ) -1001020304050 Resolution protonBias protonResolution H&HeBias H&He
Fig.
14: Bias and resolution of reconstructed energy for proton and H&He.Based on the aperture described above, the cosmic ray spectra of proton and H&He arepredicted according to the H¨orandel model [21] and the ARGO-YBJ & WFCT model [6] whichis an experimental model. The exposure time is set to one year with 10% duty cycle. Sincethe experimental model ARGO-YBJ & WFCT only shows the energy spectrum of H&He, thespectrum of pure proton is obtained by dividing by 2, that is, the ratio of proton and helium is1:1. The bending of the knee keeps unchanged.
Energy (GeV) - S r) ( m . ( G e V ) . F l u x * E ProtonH & He Energy (GeV) - S r) ( m . ( G e V ) . F l u x * E ProtonH & He
Fig.
15: The prospect of spectra at the H¨orandel model (left) and the ARGO-YBJ & WFCTmodel (right).As shown in Fig. 15, the left panel shows the spectrum with the H¨orandel model; and theright panel shows the ARGO-YBJ & WFCT model. The black dots show proton and the red dotsshow H & He . The statistical error are very small with the aperture and effective time presentedabove. The corresponding event rate in one year is shown in Fig. 16.At 800 TeV, 5121/6461 H&He events and 1130/1476 proton events can be selected in oneyear according to the H¨orandel model and the ARGO-YBJ & WFCT model respectively. At 2PeV, 1023/849 H&He events and 222/210 proton events can be selected in one year accordingto the H¨orandel model and the ARGO-YBJ & WFCT model respectively.Therefore, if the knee of the cosmic ray light component is below 1 PeV, the 1/4 LHAASOarray can give good measurements within one year. If the knee of light component is higherthan 3 PeV, more observation time or more effective methods is needed to get enough statistics.Considering the application of SiPM, it allows observation in the moon night of WFCTA. Theduty cycle of the hybrid observation can be extended. Moreover, during the construction ofLHAASO, more Cherenkov telescopes can also increase the statistics effectively . XXXXXX-14hinese Physics C Vol. XX, No. X (2019) XXXXXX
Energy (GeV) N u m be r o f E v en t s ProtonH & He
Energy (GeV) N u m be r o f E v en t s ProtonH & He
Fig.
16: The statistics of one year with 10% duty cycle at the H¨orandel model (left) and theARGO-YBJ & WFCT model (right).
The simulation for the second stage of 1/4 LHAASO hybrid detection is operated. Threetypes of detectors, WFCTA, WCDA and MD, are included in this study. Parameters from eachdetector array are studied and tuned in detail. After removing the highly correlated parameters,six parameters are used as input of the BDTG classifier for TMVA machine training and testing.The results show that the classification of pure proton is weaker than that of H&He.After cutting, high purity light component samples were selected. The aperture of pureproton is 900 m Sr , the contamination is around 10% according to the H¨orandel model; theaperture of H&He is 1800 m Sr , the contamination is less than 3%. For the aperture H&He,the uncertainty of the different strong interaction models is about ± ±
4% and the resolution is less than 20% for both pure protonand H&He samples.According to the energy spectrum expectation of the ARGO-YBJ & WFCT model, thespectrum of light component in cosmic rays will be measured accurately in a short time by 1/4LHAASO array.
This work is supported by the National Key R&D Program of China (NO. 2018YFA0404201and NO. 2018YFA0404202). This work is also supported by the Key Laboratory of ParticleAstrophysics, Institute of High Energy Physics, CAS. Projects No. Y5113D005C, No. 11563004and No. 11775248 of National Natural Science Foundation (NSFC) also provide support to thisstudy.
References γ Coll.), Astrophysical Journal, 2008, arXiv:0801.1803v1.3 M. Aglietta et al.(EAS-TOP Coll.), Astroparticle Physics, 10:l-9.4 Johannes Bl¨umer, Ralph Engel, J¨org R.H¨orandel, Progress in Particle and Nuclear Physics 63(2009) 293-338.5 T. Antoni et al., Astroparticle Physics, 2005, Vol. 24(1-2):Pages 1-25.6 B. Bartoli et al.(for ARGO-YBJ and LHAASO Coll.), Phys. Rev. D, 2015, 92(9): 25-31.7 Z. Cao (for LHAASO Coll.), Chin.Phys. C, 2010, (34)249-252.8 Z. Cao (for LHAASO Coll.), 31th ICRC, 2009.9 Peter K.F.Grieder, ISBN978-3-540-76940-8(2010).10 B. Bartoli et al.(for ARGO-YBJ Coll.), Astrop.Phys., 90 (2017) 20-27.11 L.L. Ma, 33th ICRC, 2013, p.1013.12 J.L. Liu, 33th ICRC, 2013, p.0880.13 Q. An et al.(for LHAASO Coll.), Nucl.Instrum.Meth. A, 2013, 724: 12-19.14 C. Liu et al. (for LHAASO Coll.), 35th ICRC, 2017, p.301.
XXXXXX-15hinese Physics C Vol. XX, No. X (2019) XXXXXX
15 S. S. Zhang et al.(for LHAASO Coll.), TIPP 2017.16 Baiyang Bi et al.(for LHAASO Coll.), Nucl.Instrum.Meth. A, 899 (2018) 94-100.17 S.S.Zhang et al., Nucl.Instrum.Meth. A, 2012, 629(1)57-65.18 Zhao Jing et al.(for LHAASO Coll.). Chin.Phys. C, 2014, 38(3) 036002.19 X. Zuo et al.(for LHAASO Coll.), Nucl.Instrum.Meth. A, 2015, 789: 143.20 D. Heck and T. Pierog, Extensive Air Shower Simulation with CORSIKA: A User’s Guide, 2015.21 J¨org R.H¨orandel, Astroparticle Physics, 2003, 19(2): 193-220.22 Xiurong Li et al. (for LHAASO Coll.), 35th ICRC, 2017, p.381.23 Yu.A. Fomin et al., Nuclear Physics B (Proc. Suppl.) 175-176 (2008) 334-337.24 Hillas, A. 1985, 19th ICRC Vol.3, 445-448.25 A. Hoecker, P. Speckmayer et al, TMVA Users Guide, 2013, arXiv:physics/0703039.26 P Speckmayer, A H¨ocker et al., Journal of Physics: Conference Series 219 (2010) 032057.27 B. Bi, Z. Cao et al.(for LHAASO Coll.), 2016, Frascati Physics Series Vol. 64.28 T Lamp´en, F Garcia et al., Journal of Physics: Conference Series 119 (2008) 03202829 Zizhao Zong et al.(for LHAASO Coll.), 35th ICRC, 2017(547).15 S. S. Zhang et al.(for LHAASO Coll.), TIPP 2017.16 Baiyang Bi et al.(for LHAASO Coll.), Nucl.Instrum.Meth. A, 899 (2018) 94-100.17 S.S.Zhang et al., Nucl.Instrum.Meth. A, 2012, 629(1)57-65.18 Zhao Jing et al.(for LHAASO Coll.). Chin.Phys. C, 2014, 38(3) 036002.19 X. Zuo et al.(for LHAASO Coll.), Nucl.Instrum.Meth. A, 2015, 789: 143.20 D. Heck and T. Pierog, Extensive Air Shower Simulation with CORSIKA: A User’s Guide, 2015.21 J¨org R.H¨orandel, Astroparticle Physics, 2003, 19(2): 193-220.22 Xiurong Li et al. (for LHAASO Coll.), 35th ICRC, 2017, p.381.23 Yu.A. Fomin et al., Nuclear Physics B (Proc. Suppl.) 175-176 (2008) 334-337.24 Hillas, A. 1985, 19th ICRC Vol.3, 445-448.25 A. Hoecker, P. Speckmayer et al, TMVA Users Guide, 2013, arXiv:physics/0703039.26 P Speckmayer, A H¨ocker et al., Journal of Physics: Conference Series 219 (2010) 032057.27 B. Bi, Z. Cao et al.(for LHAASO Coll.), 2016, Frascati Physics Series Vol. 64.28 T Lamp´en, F Garcia et al., Journal of Physics: Conference Series 119 (2008) 03202829 Zizhao Zong et al.(for LHAASO Coll.), 35th ICRC, 2017(547).