DsTau: Study of tau neutrino production with 400 GeV protons from the CERN-SPS

In the DsTau experiment at the CERN SPS, an independent and direct way to measure tau neutrino production following high energy proton interactions was proposed. As the main source of tau neutrinos is a decay of Ds mesons, produced in proton-nucleus interactions, the project aims at measuring a differential cross section of this reaction. The experimental method is based on a use of high resolution emulsion detectors for effective registration of events with short lived particle decays. Here we present the motivation of the study, details of the experimental technique, and the first results of the analysis of the data collected during test runs, which prove feasibility of the full scale study of the process in future.


Introduction
Tau neutrino existence was predicted after the tau lepton discovery in 1975 [1], and then it was further strengthened by the confirmation of three active neutrinos by LEP in 1989 [2][3][4][5]. However, only in 2000, 25 years later from the discover of tau lepton, first tau neutrinos were detected in the DONuT experiment [6]. In 2008, DONuT published the tau neutrino cross section result with a rather large systematic uncertainty [7]. In 2015, flavour changing neutrino oscillations were established by detecting appearance signals of neutrino oscillations, including ν τ appearance through ν µ,e ν τ oscillations by OPERA [8] and Super-Kamiokande (SK) [9]. IceCube [10] also reported an evidence of ν τ appearance recently.
Given the poor measurements of tau neutrinos, their properties are not well studied. In particular, the cross section of tau neutrino charge current (CC) interaction is known with much larger statistical and systematical uncertainties compared to the other neutrino flavors, as shown in Figure 1. On the other hand, a precise measurement of the tau neutrino interaction cross section would be interesting since it could allow testing of the Lepton Universality (LU) in neutrino scattering. LU is a principal assumption of the Standard Model (SM) of particle physics, however its validity was questioned by recent results on the B decay asymmetry [11][12][13]. There is an expectation of the same effect in ν τ scattering [14]. The measurement of the ν τ CC cross section has also a practical impact on current and future neutrino oscillation experiments. Mass hierarchy measurements in the atmospheric Super-Kamiokande (SK) [15] and accelerator neutrino experiments (for example, in DUNE [16] and HyperKamiokande [17]) rely on ν e flux measurements, but the sample of registered events is contaminated by ν τ interactions due to τ → e decays. Unlike other error sources, the uncertainty of the ν τ cross section cannot be constrained by near detector measurements. The systematic uncertainty from the ν τ interaction cross section will be a limiting factor in oscillation analyses in these experiments [18,19]. A better knowledge of the tau neutrino properties would also be a basic input to the astrophysical ν τ observation by IceCube [20].  Figure 1. Left: ν, ν averaged energy independent cross section of the three neutrino flavors ( [21] for ν e , [22] for ν µ , and for ν τ [7]). The SM prediction is indicated as a dashed horizontal line. For the DONuT result, since there is no measurement of the parameter n concerning the D s double differential production cross-section (Eq. 1.1), the value is plotted in the empirical range of a parameter n given by the DONuT paper as in the right plot. Right: The cross section result by the DONuT experiment, as a function of the parameter n [7]. The horizontal dashed line corresponds to the predicted value from lepton universality. The arrow shows the value of n estimated from PYTHIA simulation (version 6.129).
So far the tau neutrino interaction cross section was only measured in the DONuT [7], OPERA [8] and SK [23] experiments though under rather different conditions. DONuT was the first experiment which has directly registered 9 events of the ν τ CC interactions in the accelerator-based neutrino beam and extracted the value of the cross section under certain assumptions on the tau neutrinos flux. SK reported the cross section measurement with oscillated atmospheric neutrinos, based on an estimation of amount of tau neutrino events in the total sample by the Neural Network analysis. The energy ranges of the ν τ registered in both experiments are rather different, which makes the results difficult for comparison as different processes are responsible for the tau neutrino scattering. Nevertheless SK has reported on the consistency of the results with the DONuT measurements under several assumptions. OPERA registered 10 oscillated tau neutrinos and published the cross section estimation though with a rather large error [24]. In general, all the measurements performed so far have large statistical and/or systematical errors of 30-50% due to low statistics and experimental uncertainties. In a future experiment at CERN, SHiP [25], it is proposed to perform a rich neutrino program [26] and is expected to collect thousands of tau neutrino interactions, hence providing a negligible statistical error of the cross section measurement. Therefore the overall accuracy of the cross section will be determined by the systematic error, and, in particular, by the ν τ flux uncertainty, which is to be studied by the DsTau experiment.
The dominant source (> 90%) of ν τ in the accelerator-based neutrino beam is leptonic decays of D ± s mesons produced in proton-nucleus interactions: producing ν τ and ν τ in every decay. The tau neutrino beam produced in this way has rather large beam divergence due to the emission angle of D s and also the decay momenta. Taking into account the finite acceptance of the detectors, it is necessary to know the differential production cross section of D s to calculate the tau neutrino flux on the detector. For example, in the case of DONuT only about 6% of the produced ν τ intersected the DONuT detector, and 5% are expected in SHiP [25]. Conventionally, the differential production cross section of charmed particles is approximated by a phenomenological formula √ s ) and p T is the transverse momentum. n and b are the parameters controlling the longitudinal and transverse dependence of the differential production cross section, respectively. Although there were several measurements on charm particles [27][28][29][30][31], there is a lack of measurements on the D s differential production cross section in the proton interactions, especially concerning the longitudinal dependence represented by the parameter n (see a detailed review in Appendix A). This has been the main uncertainty of the ν τ cross section measurements.
At the deep inelastic interaction regime like in DONuT, the tau neutrino cross section, can be approximated as where σ const ντ is an energy independent cross section, E -the neutrino energy and K(E) -a factor taking into account kinematic effects due to the τ -lepton mass. The lack of experimental constraints on the parameter n relating to the D s differential production cross section forced the DONuT collaboration to present their result as a function of n, as σ const ντ = 2.51n 1.52 (1 ± 0.33(stat) ± 0.33(syst)) × 10 −40 cm 2 GeV −1 .
(1.3) Figure 1 on the right shows this DONuT result plotted in the empirical range of n.
DONuT also presented the numerical value for the cross section by making use of the value n = 6.1 derived from the PYTHIA simulation [32] (indicated in Figure 1), as σ const ντ = (0.39 ± 0.13(stat) ± 0.13(sys)) × 10 −38 cm 2 GeV −1 . (1.4) Note that the contribution due to the uncertainty of n is omitted from this formula. When it is included, the relative systematic uncertainty exceeds 50%.
Thus, a new measurement of differential production cross section of D s is vital for future neutrino experiments as well as the re-evaluation of the DONuT result.
In the DsTau experiment, a direct study of tau neutrino production, namely a measurement of D s → τ → X decays following high-energy proton-nucleus interactions, is proposed. DsTau exploits a simple setup consisting of a segmented high resolution vertex detector capable to recognize D s → τ → X by their very peculiar double-kink topology. The project aims to detect ∼ 1000 D s → τ → X decays in 2.3×10 8 proton interactions with the tungsten target. State-of-the-art nuclear emulsion detectors with a nanometricprecision readout will be used to achieve this goal. The modern use of the emulsion detection technology is based on the fast evolution of the high-speed and high-precision automatic readout of emulsions developed during the last two decades and available today [33][34][35].
DsTau will provide an independent ν τ flux prediction for future neutrino beams with accuracy under 10%. Then, the systematic uncertainty of the ν τ CC cross section measurement can be made sufficiently low to test LU in neutrino scattering.
In addition to the primary aim of measuring D s differential production cross section, in 2.3×10 8 proton interactions, a high yield of O(10 5 ) charmed particle pairs is expected. The analysis of those events can provide valuable by-products, such as a measurement of the intrinsic charm content in proton [36] by measuring the emission angle (pseudorapidity) of the charmed particle pairs, the interaction length of charmed hadrons, the Λ c production rate and search of super-nuclei [37], that have never measured.

Overview of project
The topology of D s → τ → X events appears as a "double-kink", as shown in Figure 2. In addition, because charm quarks are created in pairs, another decay of a charged/neutral charmed particle will be observed with a flight length of a few millimetres. Such a doublekink plus decay topology in a short distance is a very peculiar signature of this process, and the background mimicking this topology is marginal. However, to register the events is a challenge. First, all the decays take place at a scale of millimetres. Second, the kink angle of D s → τ is anticipated to be very small. The expected signal features were studied making use of Pythia 8.1 [38], as shown in Figure  3. The mean flight lengths of D s , τ and pair-charms are 3.6 mm, 2.1 mm and 4.2 mm, respectively. Although the kink angle at τ decays vertex is easily recognizable (mean kink angle of 96 mrad), the kink angle at D s → τ decays is rather small, 6.2 mrad on average. The measurement of D s momentum is also difficult as two ν τ escape measurement, spoiling the invariant mass reconstruction. Therefore we'll use an emulsion detector with nanometric spatial accuracy to efficiently recognize the events.  The DsTau module structure is shown in Figure 4. The structure is divided into two parts. The upstream part is named the decay module. The basic unit is made of a 500-µmthick tungsten plate (target) followed by 10 emulsion films interleaved with 9 200-µm-thick plastic sheets which act as a decay volume for short-lived particles as well as high-precision particle trackers. This structure (thickness of 5.3 mm) is repeated 10 times. Five additional emulsion films are placed most upstream of the module to tag the incoming beam protons. It is then followed by the downstream part made of a repeated structure of emulsion films and 1-mm-thick lead plates for momentum measurement of the daughter particles. The entire detector module is 12.5 cm wide, 10 cm high and 8.6 cm thick and consists of a total of 131 emulsion films. Momentum of the particles can be determined by Multiple Coulomb Scattering (MCS) measurement along their tracks in the tungsten plates and in the downstream ECC [39].
The emulsion detector ( Figure 5) comprises a countless number of silver halide crystals dispersed in the gelatin media. Each of those crystals is a semiconductor with a band gap of 2.5 eV and works as an individual detection channel for charged particles. The high density of detection channels, of the order of 10 14 channels/cm 3 , makes this detector suitable for three-dimensional high-precision measurements.
The emulsion films for DsTau have two emulsion layers (70 µm thick) poured onto both sides of a 210-µm-thick plastic base 1 . The diameter of the silver halide crystals is 200 nm, as shown in Figure 5 (middle). Once a charged particle passes through the emulsion layer, the ionization is recorded quasi-permanently, and then amplified and fixed by the chemical process. A trajectory of a charged particle can be observed under an optical microscope as shown in Figure 5 (right). An emulsion detector with 200-nm-diameter crystals and a 210-µm-thick base has a track position resolution of 50 nm [41] and an angular resolution . Schematic view of the module structure. A tungsten target plate, is followed by 10 emulsion films alternated by 9 plastic sheets acting as a tracker and decay volume of 4.8 mm. The sensitive layers of emulsion detectors are indicated by green color. This basic structure is repeated 10 times, and then followed by a lead-emulsion ECC structure for momentum measurement of the daughter particles. Figure 5. A cross-sectional view of an emulsion film [40] (left), electron microscope picture of silver halide crystals (middle) and a track of a relativistic particle (10 GeV/c π beam) (right). of 0.34 mrad (projection). With this angular resolution, one can detect 2-mrad kink with 4σ confidence.
A key feature of the modern emulsion technique is the use of the fast readout instruments, which allow extracting and digitizing the information on the tracks fully automatically. Emulsion detectors and the automated readout systems have been successfully employed by several neutrino experiments such as CHORUS [42], DONuT [6,7] and OPERA [8]. The latest scanning system, the HTS system [33,34], allows scanning of the emulsion films at the speed of 5,000 cm 2 per hour per emulsion layer, which is O(100) faster than those used in OPERA. The high performance is achieved by combining a custom-made wide-field objective lens and 72 high-speed image sensors. The image data is processed in real-time manner by a dedicated track recognition software based on the GPU technology. An effective scanning performance reached a 1000 m 2 /year.
The D s momentum (P Ds ) will be reconstructed by means of a machine-learning algorithm using topological variables of the event. There are two flight lengths (F L D S , F L τ ) and two kink angles (δθ Ds→τ , δθ τ →X ), those are proportional and inverse proportional to the γ factor of decaying particles. The combination of these four variables effectively provides an estimate of P Ds . The algorithm was trained and tested with the simulated sample (D s → τ → 1 charged particle) ( Figure 6). The momentum resolution is estimated to be 20%. The detection efficiency for the D s → τ → 1 prong events (85% of τ decays) is estimated by PYTHIA 8.1 [38] simulation. Table 1 shows the selection criteria and breakdown of the estimated efficiency. The overall registration efficiency of the reaction is estimated to be 20%. Figure 7 shows the estimated detection efficiency as a function of the Feynman x (x F = 2p CM Z / √ s) (left) and the x F distributions after the selection (right). Figure 8 shows the x Fp T distribution before the selection (left) and the detection efficiency in the x Fp T plane (right). The efficiency for high x F component falls by Selection (2) concerning the kink angle ∆θ Ds→τ ≥ 2 mrad, for which a further improvement of the detection efficiency is under study. Since the angular resolution depends on the track length, it is expected to be better for the particles passing through more than two sensitive layers (see also Figure  14). So, the efficiency can be further improved by applying the varying threshold for the first kink angle, depending on the flight length of the parent and daughter.

Selection
Efficiency (%) (1) Flight length of D s ≥ 2 emulsion layers 77 (2) Flight length of τ ≥ 2 emulsion layers and ∆θ Ds→τ ≥ 2 mrad 43 (3) Flight length of D s < 5 mm and flight length of τ < 5 mm 31 (4) ∆θ τ →X ≥ 15 mrad 28 (5) Pair charm: 0.1 mm ≤ flight length < 5 mm 20 (charged decays with ∆θ ≥ 15 mrad or neutral decays) Table 1. Breakdown of the efficiency estimation for D s → τ → 1-prong decay.    The main background to D s → τ → 1 prong events are hadron interactions that can mimic decays of short lived particles. The probability of background event appearance is calculated as follows, Here, P int is the probability of a primary proton to interact in the tungsten target in the decay module. P Ds→τ BG , P τ →1 prong BG and P pair charm BG are the probabilities of secondary hadron interactions to be identified as D s → τ , τ → 1 prong and pair charm decays, respectively. These probabilities are obtained by simulating 3 × 10 5 protons on the detector with the FLUKA [48,49] simulation tools. The criteria, used for charged charm or tau decay topology selection, are applied to interactions of secondary hadrons with only one charged daughter particle (P > 2 GeV/c). In addition to high energy particles a large part of interactions has associated nuclear fragments, which is strong evidence of hadron interactions. Those are effectively rejected by requesting only one charged daughter.
The sources of background to neutral pair charm appearance are neutral hadron interactions and charged particle interactions within 2 emulsion plates from the proton interaction vertex in case such short tracks are not reconstructed. These backgrounds can be associated with any proton interactions occurred upstream. To estimate the background probability to neutral charm topology, neutral particle interactions and charged particle interactions within 2 emulsion plates were selected and tested with every primary proton interaction vertex occurred upstream in tungsten. The selection criteria in Table 1 were applied to each primary-secondary interactions couple. Search for nuclear fragments associated with secondary interaction allows suppressing background in this channel as well.
The probability for an incident proton to interact in tungsten and the probabilities of background event appearance with a certain path length and kink angles θ kink and opening angle between daughters θ open used for different signal channels are shown in Table 2. The total probabilities to register background events to double kink with charged pair charm or with neutral pair charm are 1.3 ± 0.4 × 10 −9 and 2.7 ± 0.8 × 10 −9 per incident proton, respectively. The expected numbers of background events in the full statistics of DsTau (4.6 × 10 9 p.o.t.) are 6.0 ± 1.8 and 12.4 ± 3.7 for these 2 signal channels, respectively.
DsTau will provide a differential cross section of D s meson production in 400 GeV proton-nucleus interaction in 2D space, e.g. x F − P T or P − P T , so that future neutrino experiments can directly implement it. It may also be fit with the phenomenological formula, Eq. 1.1, and the parameter n is estimated, which is relevant for a re-evaluation of the tau neutrino cross section measurement by the DONuT experiment. The expected precision of the parameter n as a function of available statistics is calculated with a toy MC method and given in Figure 9. At the statistics of 1000 D s → τ → X detected events, the relative uncertainty of the ν τ flux will be reduced to below 10%.
In order to collect 1000 D s → τ events, a 230 millions of proton interactions are to be analyzed, which is another challenge from point of view of the track density and amount of data to be processed. The high proton density of 10 5 /cm 2 at the upstream surface of emulsion detector was chosen to maximize the number of interactions in a single module. The track density will then increase in the detector, yet not exceeding 10 6 /cm 2    at the downstream part of the decay module, which is affordable for the emulsion detector readout and reconstruction. With this density, 6.25×10 5 proton interactions are expected in the tungsten target in a decay module. To accumulate 2.3 × 10 8 proton interactions in the tungsten plates, 4.6 × 10 9 protons on target are needed. More than 368 modules with a total film area of 593 m 2 will be employed for this measurement.
3 Beam exposure and analysis scheme

Data taking
Two test beam campaigns were held at CERN SPS in 2016 and in 2017. In 2018, a pilot run was conducted aiming at recording 10% of the experimental data. Here we present first results of the beam tests in 2016 and 2017. The emulsion gel was made of 200-nm-diameter silver bromide crystals, whose volume occupancy were 40% at the dried state. The gel was designed by Nagoya University and a consignment production was done by FUJIFILM in Japan. 15 m 2 of emulsion films were then produced at the special facility in the University of Bern. The detector modules were assembled at CERN, right before the exposure.
A schematic view of the detector setup is shown in Figure 10, and a photo of actual setup at the H4 beamline in November 2016 is also given. The detector module was driven by a target mover so that it was uniformly exposed to the proton beam at a density of 10 5 protons/cm 2 .
The proton beam profile was measured by a silicon pixel telescope, which comprises two planes of the ATLAS IBL modules with FE-I4A front-end readout chips [44,45]. The sensor dimension is 2.04 cm × 1.87 cm, it contains 80 × 336 pixels with a size of 250 µm × 50 µm. To have uniform irradiation of the detector, the beam spot was enlarged to, for example, 1 cm× 2 cm by tuning the beam optics, as shown in Figure 11. Each emulsion detector module was mounted on a motorized X-Y stage (target mover) to change the position of the module with respect to the proton beam, so to make the detector surface uniformly irradiated at a density of 10 5 tracks/cm 2 . The target mover consists of two axes driven by stepping motors and controlled by a Raspberry PI, which receives the SPS timing signal and the scintillation counter signal. During the beam spill, the target mover moves the detector in a horizontal direction perpendicular to the beam (X direction). When the end of the detector module in X is reached, the target mover comes back to the starting point in X with an increased vertical position (Y) by 1 cm. A scintillation counter was implemented behind the target mover to monitor the rate of protons (See Figure 10). The discriminated signal was sent to a prescaler to reduce effective number of pulses, and then counted by Raspberry PI directly. The measured rate is fed into the target mover every 0.2 second and the speed in X direction was set proportional to the rate, as shown in Figure 11 on the right. Figure 12 on the right shows a comparison of the proton beam position profiles with a constant speed control during beam on and with the intensity-driven control. A better uniformity in proton irradiation was achieved with the latter one. The exposure was done at room temperature. And then the modules were immediately cooled down in the fridge (5 • C) to prevent the fading of recorded tracks. The development of films was done at the University of Bern.
The emulsion detector is both a detector and the data storage media at the same time. To perform an efficient detection of the small kinks of D s → τ decays, the analysis is performed in two stages: (1) scan the full detector by a fast system with relatively coarse angular resolution and detect events those have two decays in a short distance, namely the decays of τ and partner charm (D ± and D 0 ). (2) perform a high precision measurement around the τ decay candidates to find D s → τ small kinks. The task (1) is performed by the fast scanning system, the HTS system. The angular resolutions of the track segments read out by the HTS is about 2.5 mrad for the angular range relevant for DsTau.
The second stage of the emulsions scanning (2) is performed at dedicated scanning stations which are very precise though not as fast as HTS. In these stations, a further improvement of the resolution of the track readout in the emulsion is achieved by replacing the conventional Z-axis stage (driven by a linear rail and stepping motor) by a piezobased Z axis (PI PIFOC R P-725.4CD). This allows suppressing mechanical vibration in the horizontal plane, which limited the readout resolution in the conventional system. A reproducibility of single hit position measurement was achieved to be 8 nm. The angular measurement reproducibility is found to be 0.16 mrad (RMS), including measurement effects (such as skewing, hysteresis, encoding resolution of X/Y axes and optical aberrations).

Reconstruction and analysis
The automatic scanning systems read out the track information accumulated in the emulsion film during the exposure, digitize it and transfer to the computers for the pattern recognition and track analysis like in case of any electronic detector. The output of the readout is the information on the track segments recorded in the top and bottom layers of the film (microtracks) 2 . A segment made by linking the microtracks on two layers in a film is called a basetrack, which is a basic unit of track information from each emulsion film for later processing. Each basetrack provides a 3D coordinates X = (x, y, z), 3D vector V = (tan θ x , tan θ y , 1) and dE/dx parameter. The basic concept of track reconstruction is to link basetracks on different films by using their position and angular correspondences. The track density in DsTau (10 5 − 10 6 /cm 2 in small angular space) is relatively high. The conventional reconstruction tools for OPERA, which had a track density of 10 2 −10 3 /cm 2 in large angular space, are often not appropriate. A new tracking algorithm has been developed to overcome this problem of reconstruction in high track density in small angular space. More detail is given in Appendix B.
An example of the reconstructed data from the detector is shown in the Figure 13 (left), which shows about ten thousand tracks in 2 mm× 2 mm ( 2.5 × 10 5 /cm 2 ).
The basetrack efficiency is measured by using reconstructed tracks and shown in Figure  13 (right). The average basetrack efficiency is higher than 95%. Although it is a minor effect, there is a downward trend of the efficiency along the depth in the module. This is likely due to the increase of track density. The obtained basetrack efficiency is high enough to assure the track detection with high efficiency (>99%).
The data processing is divided into sub-volumes, for example 2 cm × 2 cm × 15-30 emulsion films. An alignment between films (rotation, transverse position shifts and gap) is obtained by means of recorded tracks. In particular for the transverse position alignment, 400-GeV proton tracks those are supposed to be most straight are selectively used. The obtained alignment precision with 15 films as a function of data processing unit diameter is shown in Figure 14 (left). The algorithm assures a sub-micron alignment in relatively large reconstruction unit. The alignment precision has a dependence on the processing area size, which is due to a non-linear distortion of the plastic base. In order to test the quality of the alignment and reconstruction, the proton beam tracks were studied in detail. The beam angular distributions measured in 2 cm × 2 cm × 15 films (7 mm thick) are shown in Figure 15 (left and middle). Approximated by Gaussian, the distributions have standard deviations of 110 µrad in XZ projection and 300 µrad in YZ. The width in angular distribution reflects the beam divergence. The beam position profile was wider in YZ projection (see Figure 11), therefore the beam divergence was also larger in YZ. In fact, Figure 15 (right) is the Y position dependence of beam angle in YZ, showing a strong pattern. Each inclined line in this plot corresponds to exposure at a given Y position of the target mover. And the positive slope of the lines (440 µrad/cm) demonstrates that the beam is diverging. From the same plot, one can see that the track angular resolution is better than 100 µrad in the large volume. When one of the lines is fitted, the standard deviation of track slopes from the fit was found to be 76 µrad. The angular resolutions with reconstructed tracks as a function of track length is shown in Figure 14 (right).
The track density in the detector module rises with the depth in the detector. The measured number of tracks is compared with the FLUKA simulation as shown in Figure  16 (left). The track density at the beginning of the decay module is about 1×10 5 /cm 2 , and increases up to 5 × 10 5 /cm 2 in the last plates. This increase is due to the proton interaction daughters and their secondary interactions, as well as the electromagnetic showers.
Note that the total amount of material in the decay module corresponds to 0.1 λ int (0.05 by tungsten target, 0.05 by emulsion and plastic films), and 1.8 X 0 (1.4 X 0 by tungsten). Although the track density increases, it has a minor effect on the reconstruction since the measurements in the emulsion detector are not just points in space but micro-vectors with precise position and direction. The reconstruction quality (mis-connections, etc) would degrade if too many tracks will coincide both in their position and angular space. However, the secondary tracks in the detector have a large variety of angles, and the track density in angular space is not in fact increasing, as shown in Figure 16 (right). The plot shows the relative angle of tracks with respect to the proton beam angle. While the tail of relatively large angles increases, the track density in the region close to the beam direction (θ = 0) decreases because of proton interactions. The reconstructed tracks are then used to find vertices. Using the tracks with the angle tanθ <= 0.4 rad, conversing pattern with 4 or more tracks is considered as a vertex. Once a vertex is found, the parental proton track which has a minimum distance to the vertex within 3 µm is searched for. If the parent proton track is found, the vertex is considered as a primary proton interaction. Figure 17 shows a distribution of Z coordinate (along the beam) of the vertexes reconstructed near by the tungsten target. An enhancement of vertices in the tungsten target is evident. One can even see the micro-structure corresponding to the emulsion layers (of higher density) and plastic bases/spacers. Figure 18 shows the measured multiplicity of charged particles at proton interactions, compared with the prediction by FLUKA. Further study will be performed on the products of proton interactions. A systematic search for the decay topology of short-lived particles is applied to the found vertices. Events with decay topology are selected by applying the criteria equivalent to those described in Table 1. The statistics of the found vertexes and events with the double decay topology observed in a sub-sample of the data are shown in Table 3. It is consistent with what is expected from the simulation for the equivalent data sample. The decay length distribution of the data and the FLUKA simulation is shown in Figure 19, which demonstrates the exponential behaviour of flight length expected for charmed particle decays.
An example of double charm event candidate found with the help of the analysis scheme is shown in Figure 20. The proton interaction occurred in tungsten, 250 µm upstream of the following first emulsion film. There are 18 charged particles at the proton interaction vertex (v 1ry ), conversing with a mean minimum distance to the vertex (impact parameter) of 1.6 µm. One of those has a kink at 3.32 mm from v 1ry . The impact parameter of kink daughter to v 1ry is 174 µm. Additionally, two charged particles start 2. 20 Table 3. Statistics found in the sub-sample of data in which 5 629 670 protons are analyzed in one unit (1 tungsten plate). The systematic uncertainty on the expected number of vertices is due to the tungsten thickness fluctuation of 3.2%. with an opening angle of 0.132 radian. The plane made of these two particles has 15.2 mrad tilted with respect to the possible neutral parent vector, meaning that this neutral decay is not two body decay. It is unlikely that the neutral decay is due to K 0 S because K 0 S mostly decays into two body.

Conclusion and outlook
The DsTau experiment is going to study tau neutrino production following high energy proton interactions, which will provide necessary information for future ν τ experiments. The letter of intent and proposal of DsTau were submitted to CERN-SPSC [46,47] in 2016 and 2017, respectively. The SPSC recommended approving DsTau in April 2019.
The test beam in 2016-2017 and pilot run in August 2018 were performed, in which we accumulated over 20 million proton interactions in the detector. The emulsion scanning and analysis of these sample are ongoing, which would allow to the re-evaluation of the ν τ cross-section by refining the ν τ flux. The first results presented here provide a convincing proof of the feasibility of the full scale study scheduled for the next physics run at CERN SPS in 2021.
6.1 ± 0.7 for the inclusive charm meson production (D ± , D 0 , D ± s ) [30]. All the other experiments using high-energy proton beams did not distinguish D s from all the other charmed particles (D ± , D 0 , Λ c ) and only provided average values of n. Because the D s produced in proton-nucleus interaction does not contain valence quarks from the initial proton, the leading particle effect is reduced. Thus, the differential production cross section should be different from other charmed particles, which could contain quarks of incoming protons. E781 (SELEX) [31] used a 600-GeV Σ − beam (instead of a proton beam) and studied the difference between D + s and D − s . The reported value for D + s , which is not affected by the leading particle effect, is n = 7.4 ± 1.0 using about 130 D + s events (within x F > 0.15) in Σ − interactions (the value for D − s , which is affected by the leading particle effect, is n = 4.1 ± 0.3). D + s in SELEX may be similar to the DONUT situation; however, the incident particle is different, the measurement is limited to x F > 0.15, and the uncertainty of n is large owing to the limited data. These results are summarized in Table 4. Results from the LHC experiments at √ s = 7, 8 or 13 TeV are not included here since the energies differ too much (400 GeV fixed target is at √ s = 27 GeV). LHCb started to measure charm productions in fixed target configuration (D 0 production was reported in [50]). LHCb might provide a good measurement also for D s in the future, however, the application of the LHCb measurements to the ν τ flux estimation is limited because the proton energy is very different (400 GeV w.r.t. 7-8 TeV); the uncertainty from the target atomic mass effect (17% relative uncertainty in DONuT [7]) is large when extrapolating from their result with noble gas targets to tungsten one. In summary, no experimental result effectively constraining the D s differential cross section at the desired level, consequently the ν τ production, exists.

B Tracking algorithm for reconstruction in high density environment
High resolution of the emulsion detector allows us to reconstruct proton interactions in an enormous pile-up of events or tracks of the order of 10 5 tracks/cm 2 in each film. This is also thank to the nature of data (basetracks, a track segments on a film), which has a vector information with the position and angular resolutions of 0.4 µm and 2 mrad, respectively. The basic concept of track reconstruction is based on the correspondence of basetracks on different films in position and angular spaces. The widely-used algorithm uses a correspondence test of two consecutive basetracks. However, it can fail in the reconstruction in high density environment like in DsTau. Especially if two or more tracks with similar direction (within a few mrad) get close to each other (in a few micron), the algorithm may not resolve the correct paths.
New tracking algorithm was developed to reconstruct tracks in the environment with high track density and narrow angular spread. When it faces multiple path candidates, it keeps all possible paths. For the case in Figure 21 Here, a average is an averaged area made by basetrack positions, a pos , and one by basetrack angle, a angle , as shown in Figure 21-(b). n is the number of basetracks involved in the path and -0.5 is an empirical value to put a higher weight on longer paths. The path with the smallest a average is chosen to be the best one. This procedure is repeated by removing the basetracks involved in the chosen paths.