text:
# On the uncertainty in the 0𝜈𝛽𝛽 decay nuclear matrix elements

## Abstract

The nuclear matrix elements $`M^{0\nu }`$ of the neutrinoless double beta decay ($`0\nu \beta \beta `$) are evaluated for 76Ge,100Mo,130Te, and 136Xe within the Renormalized Quasiparticle Random Phase Approximation (RQRPA) and the simple QRPA. Three sets of single particle level schemes are used, ranging in size from 9 to 23 orbits. When the strength of the particle-particle interaction is adjusted so that the $`2\nu \beta \beta `$ decay rate is correctly reproduced, the resulting $`M^{0\nu }`$ values become essentially independent on the size of the basis, and on the form of different realistic nucleon-nucleon potentials. Thus, one of the main reasons for variability of the calculated $`M^{0\nu }`$ within these methods is eliminated.

The discovery of oscillations of atmospheric , solar , and reactor  neutrinos shows that neutrinos have a non-vanishing rest-mass. However, since the study of oscillations provides only information on the mass-squared differences, it cannot by itself determine the absolute values of the masses (or even, at present, the sign of the mass-squared differences). Moreover, flavor oscillations are insensitive to the charge conjugation properties of the neutrinos, i.e. whether massive neutrinos are Dirac or Majorana particles. Study of the neutrinoless double beta decay ($`0\nu \beta \beta `$) is the best potential source of this crucial information on the Majorana nature of the neutrinos and on the absolute mass scale.

The search for $`0\nu \beta \beta `$ decay has a long history. So far, the decay has not been seen, but impressive limits on its half-life have been established (for the review of the field, see e.g.  or ). With the progress in oscillation studies, which now established the mass scale $`m_{scale}\sqrt{\mathrm{\Delta }m^2}`$, a wave of enthusiasm emerged in the community to develop a new generation of experiments, sensitive to such mass scale. In that context it is crucial to develop in parallel methods capable of reliably evaluating the nuclear matrix elements, and realistically asses their uncertainties. The goal of the present work is to make a contribution to that endeavor.

Clearly, the determination of the effective Majorana mass $`m_\nu `$ can be only as good as the knowledge of the nuclear matrix elements, since the half-life of the $`0\nu \beta \beta `$ decay and $`m_\nu `$ are related by

$$\frac{1}{T_{1/2}}=G^{0\nu }(E_0,Z)|M^{0\nu }|^2|m_\nu |^2,$$
(1)
where $`G^{0\nu }(E_0,Z)`$ is the precisely calculable phase-space factor, and $`M^{0\nu }`$ is the corresponding nuclear matrix element. Thus, obviously, any uncertainty in $`M^{0\nu }`$ makes the value of $`m_\nu `$ equally uncertain. In turn,

$$m_\nu =\underset{i}{\overset{N}{}}|U_{ei}|^2e^{i\alpha _i}m_i,(\mathrm{all}m_i0),$$
(2)
where $`\alpha _i`$ are the unknown Majorana phases. The elements of the mixing matrix $`|U_{ei}|^2`$ and the mass-squared differences $`\mathrm{\Delta }m^2`$ can be determined in the oscillation experiments. Using the present knowledge of these quantities, (see, e.g., ) and limits on $`T_{1/2}`$, one could decide whether the neutrino mass pattern is degenerate, or follows the inverse or normal hierarchies. If, on the other hand, the existence of the $`0\nu \beta \beta `$ decay is proven and the value of $`T_{1/2}`$ is found, a relatively narrow range of absolute neutrino mass scale can be determined, independently of the phases $`\alpha _i`$ in most situations . However, such an important insight is possible only if the nuclear matrix elements are accurately known.

The nuclear matrix element $`M^{0\nu }`$ is defined as

$$M^{0\nu }=f|\frac{M_F^{0\nu }}{g_A^2}+M_{GT}^{0\nu }+M_T^{0\nu }|i$$
(3)
where $`|i,(|f)`$ are the wave functions of the ground states of the initial (final) nuclei. The explicit forms of the operators $`M_F^{0\nu },M_{GT}^{0\nu }`$ and $`M_T^{0\nu }`$ are given in Ref. . In comparison with most of previous $`0\nu \beta \beta `$ decay studies  the higher order terms of the nucleon current are also included in the present calculation, resulting in a suppression of the nuclear matrix element by about 30 % . Note that in the numerical calculation of $`M^{0\nu }`$ here the closure approximation is avoided and the unquenched values $`g_A=1.25`$, $`g_V=1.0`$ are used.

Two basic methods are used in the evaluation of $`M^{0\nu }`$, the quasiparticle random phase approximation (QRPA) with its various modifications and the nuclear shell model (NSM). These two approaches represent in some sense opposite extremes.

In the QRPA one can include essentially unlimited set of single-particle states, but only a limited subset of configurations (iterations of the particle-hole, respectively two-quasiparticle configurations). In the context of QRPA several issues that have been raised, and deserve a systematic study:

* For realistic interactions the QRPA solutions are near the point of the so-called collapse and thus its applicability is questionable. Numerous attempts have been made to extend the method’s range of validity by partially avoiding the violation of the Pauli principle. Here we consider the simplest of them, the renormalized QRPA (RQRPA) . We compare the results of RQRPA with those of the standard QRPA .
* The choice of the size of single-particle (s.p.) space is to some extent arbitrary, often dictated by convenience. What effect this choice has on the $`M^{0\nu }`$ values is the main thrust of the present work.
* There are various forms of the nucleon-nucleon potential that lead to somewhat different forms of the resulting $`G`$-matrix ***Brueckner reaction matrix elements with these nucleon-nucleon potentials were calculated and used in Refs. .. By comparing the results obtained with three such potentials (Bonn-CD , Argonne , and Nijmegen ) we show that the resulting $`M^{0\nu }`$ are essentially identical, and independent of the choice of the realistic nucleon-nucleon potential.

In contrast, in the NSM one chooses a limited set of single-particle states in the vicinity of the Fermi level, and includes all (or most) configurations of the valence nucleons on these orbits in the evaluation of $`M^{0\nu }`$. The main open question in this approach is to determine the effects of the neglected single-particle states further away from the Fermi level. As shown below, we have also performed QRPA and RQRPA calculations with the set of s.p. states used usually in NSM. It appears that these methods, at least with the nucleon-nucleon interaction we used, are not applicable for such small s.p. bases.

In what follows the $`0\nu \beta \beta `$ decay nuclear matrix elements $`M^{0\nu }`$ for 76Ge, 100Mo, 130Te, and 136Xe are evaluated. These nuclei are most often considered as candidate sources for the next generation of the experimental search for $`0\nu \beta \beta `$ decay. For each of them three choices of the s.p. basis are considered. The smallest one has 9 levels (oscillator shells $`N`$=3,4) for 76Ge, and 13 levels (oscillator shells $`N`$=3,4 plus the $`f+h`$ orbits from $`N=5`$) for 100Mo, 130Te and 136Xe. For the intermediate size s.p. base the $`N=2`$ shell in 76Ge and 100Mo are added, and for 130Te and 136Xe also the $`p`$ orbits from $`N=5`$. Finally, the largest s.p. space  contains 21 levels for 76Ge and 100Mo (all states from shells $`N=15`$), and 23 levels for 130Te and 136Xe ($`N=15`$ and $`i`$ orbits from $`N=6`$). Thus the smallest set corresponds to $`1\mathrm{}\omega `$ particle-hole excitations, and the largest to about $`4\mathrm{}\omega `$ excitations. The s.p. energies have been calculated with the Coulomb corrected Woods-Saxon potential.

It is well known that the residual interaction is an effective interaction that depends on the size of the single-particle basis. Hence, when the basis is changed, the interaction should be modified as well. Here we propose a rather simple way to accomplish the needed renormalization.

There are three important ingredients in QRPA and RQRPA. First, the pairing interaction has to be included by solving the corresponding gap equations. Within the BCS method the strength of the pairing interaction depends on the size of the s.p. basis. As usual, we multiply the pairing part of the interaction by a factor $`g_{pair}`$ whose magnitude is adjusted, for both protons and neutrons, such that the pairing gap is correctly reproduced, separately for the initial and final nuclei.

Second, QRPA equations of motion contain a block corresponding to the particle-hole interaction, renormalized by an overall strength parameter $`g_{ph}`$. That parameter is typically adjusted by requiring that the energy of some chosen collective state, often the giant GT resonance, is correctly reproduced. We find that the calculated energy of the giant GT state is almost independent of the size of the s.p. basis and is well reproduced with $`g_{ph}1`$. Accordingly, we use $`g_{ph}=1`$ throughout, without adjustment.

Finally, QRPA equations of motion contain a block corresponding to the particle-particle interaction, renormalized by an overall strength parameter $`g_{pp}`$. (The importance of the particle-particle interaction for the $`\beta `$ strength was recognized first in Ref. , and for the $`\beta \beta `$ decay in .) It is well known that the decay rate for both modes of $`\beta \beta `$ decay depends sensitively on the value of $`g_{pp}`$. In the following we use this property to find the value of $`g_{pp}`$ for each of the possible s.p. bases. The value of the parameter $`g_{pp}`$ is fixed in each case so that the known half-life of the $`2\nu \beta \beta `$ decay is correctly reproduced. The $`2\nu \beta \beta `$ half-lives and average matrix elements collected in Table 1 of Ref.  are used, where the original references to the corresponding experiments can be found. A similar compilation of the $`2\nu \beta \beta `$ data can be found in Ref. . The resulting adjusted values of $`g_{pp}`$ are shown in Table I for both the RQRPA and QRPA methods. One can see that as the basis increases, the effective $`g_{pp}`$ decreases For 136Xe we use the lower limit on the $`2\nu `$ half-life for the adjustment of $`g_{pp}`$; hence further adjustment might be needed when the $`2\nu `$ half-life becomes known., as expected.

The adjustment of $`g_{pp}`$ is a crucial point of the present work. Several studies of the sensitivity of the nuclear matrix elements $`M^{0\nu }`$ to various modifications of the QRPA method as well as to the number of s.p. states and values of $`g_{pp}`$ were made in the recent past . Typically, these studies concluded that the values of $`M^{0\nu }`$ vary substantially depending on all of these things. Here we come to a different conclusion: by requiring that the known $`2\nu \beta \beta `$ decay half-life is correctly reproduced, and adjusting the parameter $`g_{pp}`$ accordingly, we remove much of the sensivity on the number of single-particle states, on the $`NN`$ potential employed, and even on whether RQRPA or just simple QRPA methods are used.

Clearly, the chosen procedure of finding the effective interaction is rather crude. Ideally, properly evaluated effective Hamiltonian, as well as the corresponding effective $`\beta \beta `$ operator should be used in each case. However, as shown below, the chosen procedure appears to be sufficient to stabilize the values of the nuclear matrix elements $`M^{0\nu }`$. Note, that we did not use the quenched value of the axial current coupling constant $`g_A`$, as is often done in studies of ordinary $`\beta `$ decay. We are convinced that with the quenched $`g_A`$ the basic conclusion of our work will be similar, even though the values of $`g_{pp}`$ in Table I would be of course different.

Having fixed the parameters of the effective Hamiltonian we can proceed and evaluate the $`0\nu \beta \beta `$ nuclear matrix elements $`M^{0\nu }`$ and then the corresponding half-life (we list the phase-space factors in Table II in units of $`10^{25}`$ years for $`m_\nu `$ in eV.) In evaluating $`M^{0\nu }`$ we take the short-range correlations into account in the standard way, i.e. by multiplying the operators with the square of the correlation function

$$f(r)=1e^{ar^2}(1br^2);a=1.1\mathrm{fm}^2,b=0.68\mathrm{fm}^2.$$
(4)
Note that the effect of short-range correlations reduces the matrix elements $`M^{0\nu }`$ by a factor of about two, in agreement with other evaluations.

As pointed out earlier, in the nuclear shell model an even smaller set of single-particle states is used corresponding to $`0\mathrm{}\omega `$. This choice reflects the practical computational limitations in handling the extremely large number of possible configurations, while it seems to be sufficient to describe the spectroscopy of low-lying nuclear states. In the NSM evaluation of the $`\beta \beta `$ decay rates  four s.p. orbits ($`f_{5/2},p_{3/2},p_{1/2},g_{9/2}`$) were used for 76Ge and a five orbits ($`d_{5/2},d_{3/2},s_{1/2},g_{7/2},h_{11/2}`$) for 130Te and 136Xe. These s.p. sets are free of the spurious center-of-mass states, but obviously miss a large part of the GT strength as well as of the strength corresponding to the higher multipoles. In order to describe GT transitions between low-lying states in the NSM, it is necessary to quench the corresponding strength. This is most conveniently formally achieved by using $`g_A=1.0`$ instead of the free nucleon value of $`g_A=1.25`$. We follow this prescription in our attempt to use this smallest s.p. space, and only there.

It appears that it is impossible to describe the $`2\nu \beta \beta `$ decay in such s.p. space using QRPA or RQRPA, and the nucleon-nucleon potentials employed in this work. One would have to renormalize the particle-particle block too much, with $`g_{pp}2.0`$, unlike the rather modest renormalization shown in Table I. With such large value of $`g_{pp}`$ the interaction is too far removed from the $`G`$-matrix used in the rest of this work. Therefore, one cannot expect to obtain sensible $`0\nu `$ matrix elements. In fact, we obtained very small matrix elements in this case for 130Te and 136Xe, while, perhaps accidentally, for 76Ge they are in a crude agreement with the NSM result .

We list the results with the three larger single-particle bases in Table II which represents the most significant part of the present work. As one can see by inspecting the entries, one can draw two important conclusions:

* The resulting $`M^{0\nu }`$ do not depend noticeably on the form of the nucleon-nucleon potential used. That is not an unexpected result.
* Even more importantly, with our choice of $`g_{pp}`$ the results are also essentially independent on the size of the s.p. basis. This is a much less obvious and rather pleasing conclusion. It can be contrasted with the result one would get for a constant $`g_{pp}`$ independent on the size of the s.p. basis. The values of $`M^{0\nu }`$ differ then between the small and large bases by a factor of two or more.

The effect of the $`g_{pp}`$ adjustment is illustrated in Fig. 1, showing that our procedure leads to almost constant $`M^{0\nu }`$ matrix elements. On the other hand, by choosing a fixed value of $`g_{pp}`$ the resulting $`M^{0\nu }`$ matrix elements for 9 and 21 s.p. levels would differ substantially.

The entries in Table II are relatively close to each other. To emphasize this feature, each calculated value is treated as an independent determination and for each nucleus the corresponding average $`M^{0\nu }`$ matrix elements (averaged over the three potentials and the three choices of the s.p. space) is evaluated, as well as its variance $`\sigma `$

$$\sigma ^2=\frac{1}{N1}\underset{i=1}{\overset{N}{}}(M_i^{0\nu }M^{0\nu })^2,(N=9).$$
(5)
These quantities (with the value of $`\sigma `$ in paretheses) are shown in Table III. Not only is the variance substantially less than the average value, but the results of QRPA, albeit slightly larger, are quite close to the RQRPA values. The averaged nuclear matrix elements for both methods and their variance are shown in Fig. 2.

Combining the average $`M^{0\nu }`$ with the phase-space factors listed in Table II the expected half-lives (for RQRPA and $`m_\nu `$ = 50 meV, the scale of neutrino masses suggested by oscillation experiments) are also shown in Table III. These predicted half-lives are a bit longer (particularly for the last three nuclei on our list) then various QRPA calculations usually predict. They are faster, however, then the shell model results of Ref. .

Given the average nuclear matrix elements in Table III and the phase space factors in Table II one can find a limit (or actual value) of the effective neutrino mass $`m_\nu `$ from any limit (or value) of $`T_{1/2}`$ from

$$m_\nu (\mathrm{eV})=\left(T_{1/2}\times G^{0\nu }\right)^{1/2}\times \frac{1}{M^{0\nu }}.$$
(6)
In conclusion, we have developed a “practical” way of stabilizing the values of the $`0\nu \beta \beta `$ nuclear matrix elements against their variation caused by the modification of the nucleon-nucleon interaction potential and the chosen size of the single-particle space. We have also shown that the procedure yields very similar matrix elements for the QRPA and RQRPA variants of the basic method. Even though we cannot guarrantee that this basic method is trustworthy, we have eliminated, or at least greatly reduced, the arbitrariness commonly present in the published calculations.

###### Acknowledgements.

We would like to thank Drs. H. Muether, A. Polls, and P. Czerski for providing us with the Brueckner reaction matrix elements used in the present work. V. R. was financed by the Landesforschungsschwerpunktsprogramm Baden-Wuerttemberg ”Low Energy Neutrino Physics”. The work of F. Š. was supported in part by the Deutsche Forschungsgemeinschaft (436 SLK 17/298) and by the VEGA Grant agency of the Slovac Republic under contract No. 1/0249/03. The work of P. V. was supported in part by the ”Internationales Graduiertenkolleg” GRK683/2 of the Deutsche Forschungsgemeinschaft and by the Center for Particle Physics, project No. LN00A 006 of the Ministry of Education of the Czech Republic.


format: text

text:
# MatlabMPI

## 1 Introduction

MATLAB®111MATLAB is a registered trademark of The Mathworks, Inc. is the dominant interpreted programming language for implementing numerical computations and is widely used for algorithm development, simulation, data reduction, testing and system evaluation. The popularity of Matlab is driven by the high productivity that is achieved by users because one line of Matlab code can typically replace ten lines of C or Fortran code. Many Matlab programs can benefit from faster execution on a parallel computer, but achieving this goal has been a significant challenge. There have been many previous attempts to provide an efficient mechanism for running Matlab programs on parallel computers . These efforts of have faced numerous challenges, such as limited support of Matlab functions and data structure or reliance on machine specific 3rd party parallel libraries and language extensions.

In the world of parallel computing the Message Passing Interface (MPI)  is the de facto standard for implementing programs on multiple processors. MPI defines C and Fortran language functions for doing point-to-point communication in a parallel program. MPI has proven to be an effective model for implementing parallel programs and is used by many of the worlds’ most demanding applications (weather modeling, weapons simulation, aircraft design, and signal processing simulation).

MatlabMPI consists of a set of Matlab scripts that implements a subset of MPI and allows any Matlab program to be run on a parallel computer. The key innovation of MatlabMPI is that it implements the widely used MPI “look and feel” on top of standard Matlab file I/O, resulting in a “pure” Matlab implementation that is exceedingly small ($``$250 lines of code). Thus, MatlabMPI will run on any combination of computers that Matlab supports.

The next section describes the implementation and functionality provided by MatlabMPI. Section three presents results on the bandwidth performance of the library from several parallel computers. Section four uses an image processing application to show the scaling performance that can be achieved using MatlabMPI. Section five presents our conclusions and plans for future work.

## 2 Implementation

The central innovation of MatlabMPI is its simplicity. MatlabMPI exploits Matlab’s built in file I/O capabilities, which allow any Matlab variable (matrices, arrays, structures, …) to be written and read by Matlab running on any machine, and eliminates the need to write complicated buffer packing and unpacking routines which would require $``$100,000 lines of code to implement. The approach used in MatlabMPI is illustrated in Figure 1. The sender saves a Matlab variable to a data file and when the file is complete the sender touches a lock file. The receiver continuously checks for the existence of the lock file; when it is exists the receiver reads in the data file and then deletes both the data file and the lock file. An example of a basic send and receive Matlab program is shown below

> ```
> MPI_Init;                           % Initialize MPI.
> comm = MPI_COMM_WORLD;              % Create communicator.
> comm_size = MPI_Comm_size(comm);    % Get size.
> my_rank = MPI_Comm_rank(comm);      % Get rank.
> source = 0;                         % Set source.
> dest = 1;                           % Set destination.
> tag = 1;                            % Set message tag.
> if(comm_size == 2)                  % Check size.
>   if (my_rank == source)            % If source.
>     data = 1:10;                    % Create data.
>     MPI_Send(dest,tag,comm,data);   % Send data.
>   end
>   if (my_rank == dest)              % If destination.
>     data=MPI_Recv(source,tag,comm); % Receive data.
>   end
> end
> MPI_Finalize;                       % Finalize Matlab MPI.
> exit;                               % Exit Matlab
>
> ```

The structure of the above program is very similar to MPI programs written in C or Fortran, but with the convenience of a high level language. The first part of the program sets up the MPI world; the second part differentiates the program according to rank and executes the communication; the third part closes down the MPI world and exits Matlab. If the above program were written in a matlab file called SendReceive.m, it would be executed by calling the following command from the Matlab prompt:

```
    eval( MPI_Run(’SendReceive’,2,machines) );

```

Where the machines argument can be any of the following forms:

> Run on the local host.
>
> Run on a multi processors.
>
> Run on a multi processors and communicate using dir1 and dir2, which must be visible to both machines.

The MPI\_Run command launches a Matlab script on the specified machines with output redirected to a file. If the rank=0 process is being run on the local machine, MPI\_Run returns a string containing the commands to initialize MatlabMPI, which allows MatlabMPI to be invoked deep inside of a Matlab program in a manner similar to fork-join model employed in OpenMP.

The SendReceive example illustrates the basic six MPI functions (plus MPI\_Run) that have been implemented in MatlabMPI

> Runs a matlab script in parallel.
>
> Initializes MPI.
>
> Gets the number of processors in a communicator.
>
> Gets the rank of current processor within a communicator.
>
> Sends a message to a processor (non-blocking).
>
> Receives message from a processor (blocking).
>
> Finalizes MPI.

For convenience, three additional MPI functions have also been implemented along with two utility functions that provide important MatlabMPI functionality, but are outside the MPI specification:

> Function to kill all matlab jobs started by MatlabMPI.
>
> Broadcast a message (blocking).
>
> Returns a list of all incoming messages.
>
> MatlabMPI function to prevent messages from being deleted (useful for debugging).
>
> MatlabMPI function to delete all files created by MatlabMPI.

MatlabMPI handles errors the same as Matlab, however running hundreds of copies does bring up some additional issues. If an error is encountered and the Matlab script has an “exit” statement then all the Matlab processes will die gracefully. If a Matlab job is waiting for a message that never arrives then it needs to be killed with the MPI\_Abort command. In this situation, MatlabMPI can leave files which need to be cleaned up with the MatMPI\_Delete\_all command.

On shared memory systems, MatlabMPI only requires a single Matlab license since any user is allowed to launch many Matlab sessions. On a distributed memory system, MatlabMPI requires one Matlab license per machine. Because MatlabMPI uses file I/O for communication, there must also be a directory that is visible to every machine (this is usually also required in order to install Matlab). This directory defaults to the directory that the program is launched from.

## 3 Bandwidth

The vast majority of potential Matlab applications are “embarrassingly” parallel and require minimal performance out of MatlabMPI. These applications exploit coarse grain parallelism and communicate rarely (if at all). Never-the-less, measuring the communication performance is useful for determining which applications are most suitable for MatlabMPI.

MatlabMPI has been run on several Unix platforms. It has been benchmarked and compared to the performance of C MPI on the SGI Origin2000 at Boston University. These results indicate that for large messages ($``$1 MByte) MatlabMPI is able to match the performance of C MPI (see Figure 2). For smaller messages, MatlabMPI is dominated by its latency ($``$35 milliseconds), which is significantly larger than C MPI. These results have been reproduced using a SGI Origin2000 and a Hewlett Packard workstation cluster (connected with 100 Mbit ethernet) at Ohio State University (see Figure 3).

The above bandwidth results were all obtained using two processors engaging in bi-directional sends and receives. Such a test does a good job of testing the individual links on an a multi-processor system. To more broadly test the interconnect the send receive benchmark is run on an eight node (16 cpu) Linux cluster connected with Gigabit ethernet (see Figure 4). These results are shown for one and 16 cpus. MatlabMPI is able to maintain high bandwidth even when multiple processors are communicating by allowing each processor to have its own receive directory. By cross mounting all the disks in a cluster each node only sees the traffic directed to it, which allows communication contention to be kept to a minimum.

## 4 Scaling

To further test MatlabMPI a simple image filtering application was written. This application abstracts the key computations that are used in many of our DoD sensor processing applications (e.g. wide area Synthetic Aperture Radar). The application executed repeated 2D convolutions on a large image (1024 x 128,000 $``$2 GBytes). This application was run on a large shared memory parallel computer (the SGI Origin2000 at Boston University) and achieved speedups greater than 64 on 64 processors; showing the classic super-linear speedup (due to better cache usage) that comes from breaking a very large memory problem into many smaller problems (see Figure 5).

To further test the scalability, the image processing application was run with a constant load per processor (1024 x 1024 image per processor) on a large shared/distributed memory system (the IBM SP2 at the Maui High Performance Computing Center). In this test, the application achieved a speedup of $``$300 on 304 CPUs as well achieving $``$15% of the theoretical peak (450 Gigaflops) of the system (see Figure 6).

The ultimate goal of running Matlab on parallel computers is to increase programmer productivity and decrease the large software cost of using HPC systems. Figure 7 plots the software cost (measured in Software Lines of Code or SLOCs) as a function of the maximum achieved performance (measured in units of single processor peak) for the same image filtering application implemented using several different libraries and languages (VSIPL, MPI, OpenMP, using C++, C, and Matlab ). These data show that higher level languages require fewer lines to implement the same level of functionality. Obtaining increased peak performance (i.e. exploiting more parallelism) requires more lines of code. MatlabMPI is unique in that it achieves a high peak performance using a small number of lines of code.

## 5 Conclusions and Future Work

The use of file I/O as a parallel communication mechanism is not new and is now increasingly feasible with the availability of low cost high speed disks. The extreme example of this approach are the now popular Storage Area Networks (SAN), which combine high speed routers and disks to provide server solutions. Although using file I/O increases the latency of messages it normally will not effect the bandwidth. Furthermore, the use of file I/O has several additional functional advantages which make it easy to implement large buffer sizes, recordable messages, multi-casting, and one-sided messaging. Finally, the MatlabMPI approach is readily applied to any language (e.g. IDL, Python, and Perl).

MatlabMPI demonstrates that the standard approach to writing parallel programs in C and Fortran (i.e., using MPI) is also valid in Matlab. In addition, by using Matlab file I/O, it was possible to implement MatlabMPI entirely within the Matlab environment, making it instantly portable to all computers that Matlab runs on. Most potential parallel Matlab applications are trivially parallel and don’t require high performance. Never-the-less, MatlabMPI can match C MPI performance on large messages. The simplicity and performance of MatlabMPI makes it a very reasonable choice for programmers that want to speed up their Matlab code on a parallel computer.

MatlabMPI provides the highest productivity parallel computing environment available. However, because it is a point-to-point messaging library, a significant amount code of must be added to any application in order to do basic parallel operations. In the test application presented here, the number of lines of Matlab code increased from 35 to 70. While a 70 line parallel program is extremely small, it represents a significant increase over the single processor case. Our future work will aim at creating higher level objects (e.g., distributed matrices) that will eliminate this parallel coding overhead (see Figure 8). The resulting “Parallel Matlab Toolbox” will be built on top of the MatlabMPI communication layer, and will allow a user to achieve good parallel performance without increasing the number lines of code.

## Acknowledgments

The authors would like to thank the sponsorship of the DoD High Performance Computing Modernization Office, and the staff at the supercomputing centers at Boston University, Ohio State University and the Maui High Performance Computing Center.

## Appendix A Parallel Image Filtering Test Application

```
MPI_Init;  % Initialize MPI.
comm = MPI_COMM_WORLD;  % Create communicator.
comm_size = MPI_Comm_size(comm);  % Get size.
my_rank = MPI_Comm_rank(comm);  % Get rank.
% Do a synchronized start.
starter_rank = 0;
delay = 30;  % Seconds
synch_start(comm,starter_rank,delay);
n_image_x = 2.^(10+1)*comm_size;  % Set image size (use powers of 2).
n_image_y = 2.^10;
n_point = 100;  % Number of points to put in each sub-image.
% Set filter size (use powers of 2).
n_filter_x = 2.^5;
n_filter_y = 2.^5;
n_trial = 2;  % Set the number of times to filter.
% Computer number of operations.
total_ops = 2.*n_trial*n_filter_x*n_filter_y*n_image_x*n_image_y;
if(rem(n_image_x,comm_size) ~= 0)
 disp(’ERROR: processors need to evenly divide image’);
 exit;
end
disp([’my_rank: ’,num2str(my_rank)]);  % Print rank.
left = my_rank - 1;  % Set who is source and who is destination.
if (left < 0)
  left = comm_size - 1;
end
right = my_rank + 1;
if (right >= comm_size)
  right = 0;
end
tag = 1;  % Create a unique tag id for this message.
start_time = zeros(n_trial);  % Create timing matrices.
end_time = start_time;
zero_clock = clock;  % Get a zero clock.
n_sub_image_x = n_image_x./comm_size;  % Compute sub_images for each processor.
n_sub_image_y = n_image_y;
% Create starting image and working images..
sub_image0 = rand(n_sub_image_x,n_sub_image_y).^10;
sub_image = sub_image0;
work_image = zeros(n_sub_image_x+n_filter_x,n_sub_image_y+n_filter_y);
% Create kernel.
x_shape = sin(pi.*(0:(n_filter_x-1))./(n_filter_x-1)).^2;
y_shape = sin(pi.*(0:(n_filter_y-1))./(n_filter_y-1)).^2;
kernel = x_shape.’ * y_shape;
% Create box indices.
lboxw = [1,n_filter_x/2,1,n_sub_image_y];
cboxw = [n_filter_x/2+1,n_filter_x/2+n_sub_image_x,1,n_sub_image_y];
rboxw = [n_filter_x/2+n_sub_image_x+1,n_sub_image_x+n_filter_x,1,n_sub_image_y];
lboxi = [1,n_filter_x/2,1,n_sub_image_y];
rboxi = [n_sub_image_x-n_filter_x/2+1,n_sub_image_x,1,n_sub_image_y];
start_time = etime(clock,zero_clock);  % Set start time.
% Loop over each trial.
for i_trial = 1:n_trial
  % Copy center sub_image into work_image.
  work_image(cboxw(1):cboxw(2),cboxw(3):cboxw(4)) = sub_image;
  if (comm_size > 1)
    ltag = 2.*i_trial;    % Create message tag.
    rtag = 2.*i_trial+1;
    % Send left sub-image.
    l_sub_image = sub_image(lboxi(1):lboxi(2),lboxi(3):lboxi(4));
    MPI_Send(  left, ltag, comm, l_sub_image );
    % Receive right padding.
    r_pad = MPI_Recv( right, ltag, comm );
    work_image(rboxw(1):rboxw(2),rboxw(3):rboxw(4)) = r_pad;
    % Send right sub-image.
    r_sub_image = sub_image(rboxi(1):rboxi(2),rboxi(3):rboxi(4));
    MPI_Send( right, rtag, comm, r_sub_image );
    % Receive left padding.
    l_pad = MPI_Recv( left, rtag, comm );
    work_image(lboxw(1):lboxw(2),lboxw(3):lboxw(4)) = l_pad;
  end
  work_image = conv2(work_image,kernel,’same’);   % Compute convolution.
  % Extract sub_image.
  sub_image = work_image(cboxw(1):cboxw(2),cboxw(3):cboxw(4));
end
end_time = etime(clock,zero_clock);  % Get end time for the this message.
total_time = end_time - start_time  % Print the results.
% Print compute performance.
gigaflops = total_ops / total_time / 1.e9;
disp([’GigaFlops: ’,num2str(gigaflops)]);
MPI_Finalize;  % Finalize Matlab MPI.
exit;

```


text:
# Statistical Challenges with Massive Data Sets in Particle Physics

## I Introduction

Particle physics is the science focused on the understanding of the most elementary constituents of matter and the forces governing their interactions. In stark contrast to the minute size of the objects under investigation, enormous experiments are required to further this understanding. A typical next generation experiment — the Compact Muon Solenoid (CMS)  at the Large Hadron Collider (LHC), located at the European Organization for Nuclear Research (CERN)  near Geneva, Switzerland — will weigh 12.5 thousand tonnes and have millions of channels of electronics, producing nearly 40 terabytes of data per second. A cartoon of the CMS detector is provided in Fig. 1. These data will be analyzed in real time and reduced to approximately ten terabytes per day that will be stored for later analysis. Of the one billion particle collisions occurring each second at the LHC, only a few are of interest. Finding these interesting — but possibly unanticipated — collisions in such a massive data stream represents a challenging test of forefront technology and computational algorithms.

Particle physics data come in discrete packets referred to as events. In the case of a proton collider experiment, an event is a head-on collision of two protons, each moving at nearly the speed of light. The energy of the two colliding protons can allow the creation of new matter, by Einstein’s $`E=mc^2`$, which is spewed forth as debris from the collision. Particle physics detectors are designed to record and identify this debris.

Sections II and III describe the real-time reduction of the massive data stream emerging from particle physics detectors to a handful of the most important components of those debris, each carrying a direction and energy. The total number of observables in the final state after this massive reconstruction is roughly two dozen. A cartoon characterization of these observables is presented in Fig. 2.

Sections IV and V describe two methods that use the resulting two dozen high-level observables to connect to the underlying physical theory. Our present understanding of that underlying model of particle physics is encapsulated in the quantum mechanical “Standard Model,” which predicts (probabilistically) the debris that will emerge. In practice the calculations are typically done using Feynman diagrams, named after the late Caltech professor, which serve as both a convenient calculational tool and an intuitive graphical aid. An example of such a diagram is shown in Fig. 3. A set of Feynman rules associates each piece of the diagram — each vertex and each line — with an algebraic expression, which after some manipulation allows a prediction to be made for the fraction of collisions that will produce certain types of debris.

The statistical analysis of these data occurs at many levels and in various forms. This paper sketches a few of the challenges currently faced in particle physics for which statistical methods are being employed, emphasizing examples from pattern recognition and signal identification. The ultimate challenge is to cull interesting and unanticipated events from billions of collisions in a data stream flowing at terabytes per second.

## II Triggering

If the detector components that record the debris of each collision are thought of as the sense organs of the detector, the trigger should be thought of as a hardwired reflex. Many animals are programmed to respond to quick motion in their peripheral vision; analogously, today’s collider detectors are designed to respond to events with certain characteristics. While motion in an animal’s periphery may take hundreds of milliseconds to trigger the appropriate muscular response for identifying or fleeing from an apparent threat, the experiments at the LHC must decide 40 million times per second whether the collisions that have occurred are sufficiently interesting to be worth recording permanently.

The large volume of data produced in today’s particle physics detectors precludes a full analysis of all events; it is not even physically possible to record so much data. A fast analysis of the data is therefore performed in real time at various levels of detail, and data deemed sufficiently interesting are saved for a more thorough subsequent analysis.

### II.1 Trigger levels

Many of today’s particle physics experiments employ a trigger system with three levels:

* Level 1. Custom electronics located physically on the detector perform a fast analysis within the detector components themselves, reducing the data that is subsequently passed through communication networks downstream.
* Level 2. Detector components are read out into computers that perform a regional analysis of the data.
* Level 3. Data from different regions of the detector are brought together in a single CPU to allow a final global decision.

An animal’s response to a perceived threat is determined by a neural system that can be placed somewhere between “hardware” (built-in, unmodifiable features that come with the biological package) and “software” (modifiable logic learned as the organism develops). Level 1 of many particle physics triggers make heavy use of “firmware,” often in the form of field programmable gate arrays (FPGAs). These devices combine the hardware advantage of high operational speeds with the flexibility and reprogrammability of software, allowing physicists to change the algorithms used as experience is gained.

The Level 2 trigger uses the longer decision time available to it to combine information from several sub-detectors, and to perform a more sophisticated selection of the data than is possible at Level 1. In cases for which two bytes of information are available from a particular detector element, the Level 1 trigger may use only the coarse information available in the first byte; at Level 2 the full information may be used. Calibration corrections that are not available in Level 1 can be applied to the data at Level 2. It is also possible at this stage to consider multivariate analysis methods.

At Level 3, the entire event is available for consideration. At this stage a fast version of the final analysis can be performed, and time-intensive algorithms can be brought to bear. Level 3 typically is allocated $`1`$ second of decision time per event.

### II.2 Example

The end cap muon trigger for the CMS experiment  is an example of a practical trigger architecture. This trigger identifies muon particles traveling through the end caps of the detector , the signature for which is a charged particle traveling through the entire apparatus. In the end cap region there are chambers, labeled CSC’s in Fig. 4, which can capture the ionization trail of the particle.

As can be seen in Fig. 5, a muon takes a complicated path through the chambers, due to the varying magnetic field in the steel structure in which the chambers are embedded. The Level 1 trigger’s goal is to select muons that follow paths through the detector that correspond to high momentum. Low momentum muons, which exhibit more bending as they go through the chambers, are of lesser interest. There is a probability distribution of possible bend angles and paths through the detector for a muon of given momentum. The trigger uses these distributions, generated by simulation, to select events that are likely to contain high momentum muons. The two key enabling technologies for this effort have general applicability to the real time analysis of massive streams of data: sorting in field programmable gate arrays, and the use of memory lookups.

Field programmable gate arrays are large scale chips, comprising many gates and memory cells, that can be configured to carry out custom algorithms. Gate arrays with a million logic units are not uncommon these days, and a custom CPU can be built by appropriately configuring the gates. In designing high-speed analyses of massive data streams, algorithms involving loops are to be avoided at all costs. Using FPGAs implementing massively parallel algorithms, a printed circuit board has been built that can sort in real time the nearly 30 gigabytes of data per second incident to it as part of the CMS end cap muon trigger . Such techniques have become commonplace in experimental particle physics, and have been used in a variety of triggers.

A second key enabling technology involves memory lookup devices. In the CMS end cap muon system, the muon’s momentum between two chambers can be found as a function of the difference in bend angle in those chambers,

$$p_T=\frac{A_{ij}}{\mathrm{\Delta }\varphi _{ij}},$$

where $`A_{ij}`$ is determined by fitting the results of a simulation, and allowed to vary with position in the detector. Maximization of the likelihood

$$=\frac{1}{2\pi \sigma _{12}\sigma _{23}\sqrt{1\rho ^2}}\mathrm{exp}\left(\frac{1}{2(1\rho ^2)}\left[\frac{(\mathrm{\Delta }\varphi _{12}\mu _{12})^2}{\sigma _{12}^2}\frac{2\rho (\mathrm{\Delta }\varphi _{12}\mu _{12})(\mathrm{\Delta }\varphi _{23}\mu _{23})}{\sigma _{12}\sigma _{23}}+\frac{(\mathrm{\Delta }\varphi _{23}\mu _{23})^2}{\sigma _{23}^2}\right]\right)$$
(1)
as a function of $`p_T`$ is desired, where $`\mu _{ij}=A_{ij}/p_T`$, $`\sigma _{ij}`$ is the error in $`\mathrm{\Delta }\varphi _{ij}`$, and it is assumed that the ratios $`\sigma _{ij}/\mu _{ij}`$ are independent of $`p_T`$. The result is shown in Fig. 6. The mapping between measured angle differences $`\mathrm{\Delta }\varphi _{12}`$ and $`\mathrm{\Delta }\varphi _{23}`$ and the most likely momentum $`p_T`$ is then loaded into a memory chip. This enables the calculation of a fairly complex function in the time it takes to access the memory, in this case about 12.5 nanoseconds . The computational and graphical statistics community could make welcome contributions to the determination of suitable functions for memory lookups in similar problems.

## III Reconstruction

For each one-in-a-million collision selected by the trigger, the primary components of the debris spewed forth in the collision are reconstructed. This process is analogous to the investigation of a head-on automobile collision, noting the resultant trajectories of broken doors, shattered glass, and strewn metal. The two examples considered here are track finding and jet clustering.

### III.1 Track Finding

Many of the particles produced in an event are electrically charged. These can be detected by any of several technologies, including the trail of ionization left in a volume of gas, the electron-hole pairs liberated in a semiconductor, and the scintillation light emitted upon excitation of certain complex organic molecules. In all cases the detector is trying to measure in three dimensions the passage of charged particles, usually by combining three two-dimensional views at three different angles through the detector. Tracking systems are typically embedded in a magnetic field, which bends the track by an amount inversely proportional to the momentum of the particle perpendicular to the field, allowing a measurement of this momentum.

In Fig. 7, the small dots represent the measured points of passage of charged particles. The superimposed lines are fits to the tracks. There are two stages to making these fits: pattern recognition, in which points that represent a common track are linked together; and fitting, in which track parameters, such as the momentum of the particle and its point of closest approach to the collision point, are extracted. The process is complicated by the presence of spurious noise hits, measurement error, and possible change of the particle’s direction through interaction with material. The determination of charged particle trajectories can be the most time consuming part of the reconstruction of an event; new, more efficient methods of track finding and fitting are continuously being sought.

A typical approach to track fitting involves finding a set of hit points that represent a possible track, forming a $`\chi ^2`$ quantifying how well the measured points fit particular track parameters, and then minimizing this $`\chi ^2`$ with respect to those parameters. Correlations among measurements, introduced for example by the interaction of the particle with the material it is traversing, complicates and slows the inversion of the covariance matrix required for this minimization. The Kalman filter  is commonly used in practice.

### III.2 Jet Clustering

A second common reconstruction problem in particle physics is jet clustering in calorimeters, detectors that measure the energy of particles entering them. An energetic, strongly-interacting particle entering a calorimeter causes a chain reaction of nuclear breakup and particle production, resulting in a shower of particles passing through the detector. Figure 8 shows a calorimeter unfolded out into a two-dimensional array; the height of each cell is proportional to the amount of energy deposited in it in this particular event. Quarks produced in the original collision fragment into many particles, which appear in the detector as a collimated “jet” of energy flow. The pattern recognition problem involves identifying the clumps of energy in the calorimeter that constitute a jet, allowing an estimate of the energy and direction of the parent quark. Less-than-perfect theoretical understanding of the fragmentation process complicates this already messy problem.

There are other sources of jets that one would like to identify. For example, $`\tau `$ leptons can decay into a tightly collimated jet of energy in the calorimeter. While the statistical distribution of energy from a $`\tau `$ decay is different on average from that obtained from the fragmentation of quarks, there is a significant overlap in the distributions. Optimal use of measured information to best separate jets from quarks and jets from $`\tau `$ leptons is an outstanding problem in proton collider experiments.

## IV Kernel Density Estimation

Section II discussed how potentially interesting events are triggered on in real time; Sec. III sketched how the millions of channels of detector information characterizing each event are reduced to roughly two dozen numbers describing the trajectories of the primary components of the debris from each collision. This section and the next summarize two methods designed to facilitate the interpretation of these events in the context of an underlying physical theory.

The problem of interpretation can often be formulated as a classification or discrimination problem, with the goal being to distinguish between two different hypotheses. Our understanding of the underlying physics is such that one of the two hypotheses is invariably the Standard Model, and the second is often an extension of the Standard Model. This leads naturally into thinking of the problem as searching for a “signal” over a “background.” Here and below the signal will be denoted by $`S`$, and the background will be denoted by $`B`$. The two hypotheses being compared are then schematically $`A=S+B`$ versus $`B`$. Particle physicists often perform this comparison in a rather simplistic way, carving out one contiguous region in a low-dimensional variable space, counting the number of events seen in the data within that region, and comparing to the number of events predicted from the two hypotheses $`A`$ and $`B`$ under consideration.

The modeling of the response of our detector to this debris is done through a rather detailed Monte Carlo simulation of each detector subcomponent. In the simulation the debris is tracked through detector components, causing signals on individual channels in the same way that real debris produces such signals in our physical detectors. The simulated events are then processed using the same algorithm as was used on the data. The simulation of these events in the detector often requires substantial time, on the order of several seconds per event on a modern (at the time of this writing) 2 GHz CPU. This limits the number of events available to model the predictions from different hypotheses of the underlying physics. In nearly all cases, time and other practical constraints limit the number of simulated events with which we can populate the space of two dozen observables to on the order of one million.

### IV.1 The problem

The issue described in this section could be avoided if we had an analytic form for the theoretical prediction in our space of observables. The Monte Carlo nature of our simulation prevents this; some prescription for smoothing out the density of points obtained from this Monte Carlo approach is therefore required. Even one million events do not adequately populate a ten-dimensional space, so physicists are forced to consider a subset of relevant variables for testing hypotheses of interest. The nature of each hypothesis to be tested typically immediately suggests a variable space $`𝒱`$ of low dimensionality that is apt to be most useful for distinguishing the hypothesis from the Standard Model, which has been the default hypothesis in the field for the past two decades. The dimensionality of $`𝒱`$ is typically between one and four, each dimension representing a simple algebraic combination of the original two dozen observables. In practice, the choice of which variables to use to test a given hypothesis is often a matter of much consternation and discussion.

### IV.2 Kernel solution

There are a number of ways to determine the region in which to do this counting. One such procedure, still used to a surprising extent, is to eyeball the predicted distributions for the two hypotheses and choose some rectangular box within the variable space by considering one-dimensional projections. A more sophisticated option, and one enjoying increasing popularity among particle physicists, is to use a neural network trained to distinguish between the two different hypotheses.

A third option, and the one we advocate, is to use kernels to estimate probability densities for the two hypotheses within the variable space $`𝒱`$. The simplest estimate gives the probability densities $`p(\stackrel{}{x}|S)`$ and $`p(\stackrel{}{x}|B)`$ at a point $`\stackrel{}{x}`$ within $`𝒱`$ in terms of the Monte Carlo points $`\stackrel{}{x}_{Si}`$ and $`\stackrel{}{x}_{Bj}`$ in $`𝒱`$, each with weight $`w_{Si}`$ and $`w_{Bj}`$, drawn from the probabilistic predictions of $`S`$ and $`B`$, as

$$p(\stackrel{}{x}|S)=\frac{1}{_iw_{Si}}\underset{i}{}\frac{w_{Si}}{\sqrt{(2\pi )^d\left|\mathrm{\Sigma }_S\right|}h_{S}^{}{}_{}{}^{d}}e^{\left((\stackrel{}{x}\stackrel{}{x}_{Si})^T\mathrm{\Sigma }_S^1(\stackrel{}{x}\stackrel{}{x}_{Si})/2h_{S}^{}{}_{}{}^{2d}\right)}$$
(2)
and

$$p(\stackrel{}{x}|B)=\frac{1}{_jw_{Bj}}\underset{j}{}\frac{w_{Bj}}{\sqrt{(2\pi )^d\left|\mathrm{\Sigma }_B\right|}h_{B}^{}{}_{}{}^{d}}e^{\left((\stackrel{}{x}\stackrel{}{x}_{Bj})^T\mathrm{\Sigma }_B^1(\stackrel{}{x}\stackrel{}{x}_{Bj})/2h_{B}^{}{}_{}{}^{2d}\right)}.$$
(3)
Here $`\stackrel{}{x}_{Si}`$ denotes the position within the $`d`$-dimensional variable space $`𝒱`$ of the $`i^{\text{th}}`$ Monte Carlo event used to model hypothesis $`S`$, and $`\stackrel{}{x}_{Bj}`$ denotes the position of the $`j^{\text{th}}`$ Monte Carlo event used to model hypothesis $`B`$. The width of the Gaussian kernels are set by the covariance matrices $`\mathrm{\Sigma }_S`$ and $`\mathrm{\Sigma }_B`$ of the points $`\{\stackrel{}{x}_{Si}\}`$ and $`\{\stackrel{}{x}_{Bj}\}`$, and by smoothing parameters $`h_S`$ and $`h_B`$ that depend upon the number of Monte Carlo events available and the dimensionality of $`𝒱`$.

The two densities $`p(\stackrel{}{x}|S)`$ and $`p(\stackrel{}{x}|B)`$ are then used to define a discriminant

$$D(\stackrel{}{x})=\frac{p(x|S)}{p(x|S)+p(x|B)}.$$
(4)
The discriminant $`D(\stackrel{}{x})`$ is a function of position within the variable space $`𝒱`$; at each point within that space $`D(\stackrel{}{x})`$ is a number between zero and unity, taking on values close to zero in regions that are predominantly populated by the types of collisions predicted by the hypothesis $`B`$, and taking on values that are close to unity in regions of $`𝒱`$ that are predominantly populated by the types of collisions predicted by the hypothesis $`S`$.

The set of all points $`\stackrel{}{x}`$ within $`𝒱`$ for which $`D(\stackrel{}{x})>D_{\text{cut}}`$ thus defines a region in the variable space $`𝒱`$ in which the contribution of collisions from the signal hypothesis $`S`$ is enhanced relative to the Standard Model background $`B`$. Here $`D_{\text{cut}}`$ is some number between zero and one, chosen to optimize a figure of merit. A common choice for this figure of merit is the number of events predicted by $`S`$ in the region divided by the square root of the number of events predicted by $`B`$ in the region, although more sophisticated choices for this figure of merit have been considered.

The number of events observed in the data within the region of $`𝒱`$ thus chosen is compared to the number of events predicted within that region by the “background” hypothesis $`B`$ and the “signal + background” hypothesis $`S+B`$, and from this comparison constraints are set on parameters internal to the hypothesis $`S`$, in some cases ruling out $`S`$ entirely.

One problem shared among all of these methods is that they are not sufficiently prescriptive. The simple selection of a rectangular box in a multivariate space typically takes some fiddling, since a completely prescriptive procedure for choosing this box is lacking. The parameters of a neural network typically have to be manipulated to achieve reasonable performance. In the use of kernels, the choice of the width of each kernel and its covariance also lacks a sufficiently robust and general prescription, despite many proposals existing in the literature. Obtaining a prescription that works reasonably and without pathologies in all cases would be incredibly helpful.

### IV.3 Example: Quaero

An attempt at making the procedure sketched above completely prescriptive has been implemented in an algorithm called Quaero . Particle physics experimenters are notoriously guarded with their data, rarely sharing data in “raw” form with physicists outside their large collaborations. The rationale for this is that the data in raw form are complicated by detailed effects and foibles of the detectors. In order to perform a meaningful analysis of the data, the detector must be understood in much greater detail than anyone outside of the collaboration can reasonably hope to achieve.

This presents an obstacle to progress in particle physics, since the normal steps of testing any particular hypothesis within the collaboration (and then obtaining official collaboration approval of the result) typically takes greater than a year. The solution to this conundrum is to package the knowledge of the physicists within the collaboration into a tool that can perform an analysis automatically — a tool that “knows” how to correct for the foibles of the experimental apparatus.

Figure 9 shows the Quaero web page, through which the DØ collaboration has made a subset of its data publicly available. A physicist can provide a file describing the type of events his model predicts should be seen in the detector, together with at most three variables that can be used to distinguish his signal ($`S`$) from the Standard Model ($`B`$). Quaero then constructs kernel density estimates $`p(\stackrel{}{x}|S)`$ and $`p(\stackrel{}{x}|B)`$ of the signal and background in the variables that the user has provided, uses these to construct the discriminant $`D(\stackrel{}{x})`$ according to Eq. 4, chooses a number $`D_{\text{cut}}`$ to maximize a reasonable figure of merit, and from this determines the region of the variable space $`𝒱`$ for which $`D(\stackrel{}{x})>D_{\text{cut}}`$. The number of events observed within that region in the data are then counted and compared to the number of events predicted by the physicist’s signal and to the number of events predicted by the Standard Model alone. From this comparison, parameters internal to the physicist’s model are constrained. Examples of the densities formed in this procedure are shown in Fig. 10.

Quaero uses kernels whose shape is determined from the computation of a local covariance matrix, which appears to provide reasonable performance in most of the analyses that have been attempted. It should be noted, however, that the reliance upon the “correctness” of the density estimate is in some sense minimal — the density estimate is used only within a prescription to determine a particular region within the variable space $`𝒱`$ that enhances the contribution of signal relative to background. Using the kernel density estimate directly to compute a likelihood requires confidence that the density estimate prescription used will yield non-pathologic results in all cases, and a prescription generating sufficient confidence is currently lacking.

By automating the testing of individual hypotheses, Quaero dramatically decreases the time required to test models against large sets of data. Each such test, formerly the subject of a Ph.D. theses, can now be done in an hour.

## V Sleuth

Because Quaero is designed to facilitate the testing of specific hypotheses against a large set of data, it requires that specific hypotheses be defined. A separate algorithm is required to allow searches in the data for signatures of a more general type. Sleuth looks for anomalous deviations in the data that could indicate the presence of interesting physics, allowing the simultaneous testing of a large number of vaguely-specified hypotheses .

### V.1 The problem

The Standard Model of particle physics has so far passed the experimental tests to which it has been subjected. Nonetheless, there is a hole in the Standard Model that indicates there are likely to be new fundamental discoveries at energy scales that our accelerators are just beginning to probe. Exactly what form that new physics might take is currently unknown, despite the work of hundreds of people over the last two decades.

Possible findings include

* particles with bare magnetic charge (“magnetic monopoles”);
* a new symmetry of Nature (“supersymmetry”);
* a new strong force (“technicolor”);
* the presence of a new weak force (“heavy gauge bosons”);
* the existence of large extra spatial dimensions, curled up on scales smaller than 1 mm (“extra dimensions”);
* electrons in excited states (“excited fermions”);
* laws of nature that are not isotropic (“non-commutative theories”);
* a heavier analog of the electron (“fourth generation of fermions”); and
* evidence pointing to a unification of the strong, weak, and electromagnetic forces (“grand unified theories”).

Each of these possibilities represents a class of theories with a number of adjustable parameters. When taken together, the model space is sufficiently large that systematically checking all possibilities is not a viable option. We know only vaguely what it is we should be searching for; equivalently, we are searching for more things than can possibly be tested at one time.

Humans are notoriously good at finding patterns, particularly when dealing with small numbers of events — the history of particle physics is strewn with cases of patterns being mistakenly discerned. An algorithm is required that will enable a general search for all possibilities simultaneously, rigorously taking into account the “trials factor” (a measure of the number of different possibilities considered) when reporting a final number.

### V.2 Solution

The solution of this problem in the context of particle physics at accelerators that collide protons or their antiparticles is an algorithm called Sleuth. Consideration of the many possibilities just mentioned naturally leads one to ask whether there is any common feature among them. If such a common feature exists, perhaps it can be searched for in a general way.

It turns out that such a commonality does in fact exist, justifying the following three assumptions:

* the data can be partitioned in such a way that a signal will predominantly appear within a single category;
* interactions signaling the presence of new physics will generally produce objects with larger energy than expected from background processes; and
* a signal is apt to appear as an excess of events — i.e., more events observed in the data than expected from background.

Sleuth begins by taking all of the data collected in the experiment and partitioning it into categories. This partitioning is orthogonal; each event ends up in one and only one of these categories. The partitioning is chosen such that if new physics appears in the data, it is likely to end up predominantly in a single category. Each category contains a set of “similar” events, in the sense that the events in each category contain the same types of debris from the collision.

Within each category, $`d`$ variables (where $`d`$ ranges between one and four, depending upon the category) are chosen. In the case of Sleuth, these variables correspond roughly to the energies of the objects produced in the collision. The variables to be used for each category are set by a rule agreed upon before the data are analyzed.

An arbitrary number of regions can be drawn in the variable space just defined. In order to discretize the regions that can be considered, the following procedure is adopted. A multivariate change of variables is made, such that in the new variables the prediction of the default hypothesis (the Standard Model) is a uniform distribution in $`[0,1]^d`$. We refer to $`[0,1]^d`$ for arbitrary $`d`$ as the “unit box,” meaning unit interval when $`d=1`$, unit square when $`d=2`$, unit cube when $`d=3`$, and so forth. Regions are then defined about sets of data points using the concept of Voronoi diagrams, an example of which is shown in Fig. 11. The region around any single data point is the space within the unit box closer to that data point than to any other data point in the unit box. The number of regions — all possible unions of these individual regions — that can be considered has been reduced by this process from uncountably infinite to a mere $`2^{N_{\text{data}}}`$, where $`N_{\text{data}}`$ is the number of events observed in the data in this category. Because not all such regions are physically interesting, the regions considered are restricted to those containing the “upper right-hand corner” ($`1^d`$) of the unit box and satisfying other criteria favoring regions in which all variables take on values close to 1. New physics is expected to appear in collisions with energetic particles in the final state; our mapping into the unit box is done such that these events will take on values close to the upper right hand corner.

The “interestingness” of each region is quantified as the Poisson probability that the background within that region will fluctuate up to or beyond the observed number of events. A search heuristic is used to find the most interesting region $``$ in each category in the data. The fraction $`𝒫`$ of hypothetical similar experiments, in which a set of “pseudo” data points are drawn from the background distribution, that give rise to a result more interesting than the most interesting region observed in the data is determined. $`𝒫`$ is a number between zero and unity; $`𝒫`$ will be a (uniformly-distributed) random number between zero and unity if the data are composed of background only; $`𝒫`$ should be close to zero if there really is something interesting in the data. One value of $`𝒫`$ is determined for each category. The minimal $`𝒫`$ in the (roughly 100) categories considered determines $`\stackrel{~}{𝒫}`$, the fraction of hypothetical similar experiments in which a more interesting region would be found in any of the categories considered. This $`\stackrel{~}{𝒫}`$ is the final measure of “interestingness.” If the data are distributed according to background only, $`\stackrel{~}{𝒫}`$ will be a number randomly distributed in the unit interval. If the data contain a hint of a signal, $`\stackrel{~}{𝒫}`$ will (hopefully) be small. The output of Sleuth is thus twofold: the most interesting region $``$ in any of the categories considered, and a quantitative measure $`\stackrel{~}{𝒫}`$ of how interesting that region actually is.

### V.3 Results

By construction, Sleuth will not find anything if there is nothing there to be found. But would Sleuth find something if there were something there to be found? This is impossible to address adequately in the general case; the answer depends to what degree the underlying physical assumptions made in Sleuth are realized in Nature. The question can, however, be definitively answered in any specific scenario, and several studies have been performed in order to test Sleuth’s sensitivity to a variety of new physics. The results of one such example are shown in Fig. 12. Other studies in various categories with different signals indicate that Sleuth performs as anticipated, the success of the algorithm in each case depending upon the strength of the signal and the extent to which the signal satisfies the three assumptions on which the algorithm is based. Such studies indicate that Sleuth has a reasonable shot of discovering new physics were it present in our data.

Sleuth has been applied to a range of particle collider data collected between 1992–1996 by the DØ experiment at the Fermilab Tevatron. Disappointingly, a null result has been obtained. The $`𝒫`$’s obtained in thirty-two different categories are histogrammed in Fig. 14(b), and the conversion to $`\stackrel{~}{𝒫}`$ given the number of categories considered is shown in Fig. 14(a).

### V.4 Remarks

The general problem of looking for a signal that is only vaguely specified is of course not peculiar to particle physics, and the general thread of the solution should therefore be applicable to any problem of similar type. The individual steps of the algorithm — partitioning the data into categories, defining a small set of variables for each category, rigorously defining regions within each variable space, specifying criteria for the regions to be considered, quantifying the interestingness of any particular region, searching the data to find the most interesting region within each category, conducting hypothetical similar experiments to determine $`𝒫`$ for each category, and conducting a second set of hypothetical similar experiments to determine $`\stackrel{~}{𝒫}`$ for the data as a whole — are generally modifiable to many problems. Where the size of the data set is an issue, a specific solution can often be found that scales with the number of events $`n`$ as $`𝒪(n\mathrm{log}n)`$.

## VI Summary and Discussion

The massive data sets from this decade’s particle physics experiments present a variety of challenges. The selection of which collisions to record on permanent storage (triggering), discussed in Sec. II, presents a need for smart algorithms whose implementation in a parallel hardware architecture can take the high rate of incident data. Understanding and optimizing the necessary trade-offs imposed by the constraints of each system involve design patterns that are still being learned. Imposing deadlines and ignorance of theoretical best practice result in systems that can be noticeably suboptimal.

The reconstruction of the trajectories of the primary components of the debris emerging from the collision, sketched in Sec. III, involves a number of pattern recognition problems that can be handled at the more leisurely pace of one second per collision. Track finding, a connect-the-dots problem of identifying curved particle trajectories in the presence of noise and complicating interactions of the particle with the traversed material, is prototypical problem. Tracking flying objects with sweeping radar or swimming objects emitting infrequent signals are two of many examples where similar issues are encountered. Jet clustering, the clustering of the debris observed in the detector, has similarities to many of the clustering problems currently treated in the statistics literature. These are two areas in which the statistical community’s impact on current practice in our field has been significantly less than it could be.

Kernel Density Estimation has been recognized by the field as a useful technique for estimating parent densities of finite samples, the limited statistics of which are often imposed by the time required to simulate what would be observed in our detectors if a particular model were realized in nature. Quaero, the intelligent web interface to DØ proton anti-proton collision data described in Section IV, serves as an existence proof that the testing of particular hypotheses against high energy collider data can be formalized and automated. Outstanding issues here include whether a completely prescriptive kernel density algorithm can be achieved that avoids pathologies when presented with outliers, delta functions in an otherwise continuous distribution, and discontinuous boundaries. To be a bit more provocative: For every existing kernel estimation prescription of which we are aware, a sample exists for which the prescription produces a clearly undesirable estimate of the parent density. Presumably a prescription that is never obviously unreasonable can be achieved, and this would be of great value. Such a prescription whose time cost grows only linearly with the size of the sample (even in the multidimensional case) would be of even greater value.

The problem of how to look for a vaguely specified feature in a large set of data, discussed in Section V, is omnipresent in scientific research. Having the ability to rigorously quantify the “interestingness” of a particular subset of collisions in an unbiased manner is found to be crucial. This basic lesson can be widely generalized. The solution (Sleuth) is obtained by constructing a rigorous quantification of “interestingness” before the data is analyzed, so that pseudo experiments can be simulated to determine what fraction of them would yield data as interesting as what is actually observed. Obvious applications range from identifying statistically significant galaxy clusters to identifying statistically significant disease clusters.

Opportunities for the computational and graphical statistics community to contribute to the solutions of problems faced in particle physics abound; this article has barely scratched the surface. We look forward to a productive exchange of ideas over the coming decade.

###### Acknowledgements.

P. P. thanks David Scott for the invitation to a lively conference at which this talk was given, and his colleagues at CMS and DØ for the collaborative work from which the discussion of triggers in Sec. II has been drawn. B. K. thanks his colleagues in the DØ collaboration, discussions with whom led to the development of Quaero and Sleuth, and his collaborators on CDF. P. P. is supported by a grant from the Department of Energy; B. K. is supported by a Fermi/McCormick Postdoctoral Fellowship at the University of Chicago.