MEASUREMENT OF ISOBARIC ANALOGUE RESONANCES OF 47 Ar
WITH THE ACTIVE-TARGET TIME PROJECTION CHAMBER
By
Joshua William Bradt

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics – Doctor of Philosophy
2017

ABSTRACT
MEASUREMENT OF ISOBARIC ANALOGUE RESONANCES OF 47 Ar
WITH THE ACTIVE-TARGET TIME PROJECTION CHAMBER
By
Joshua William Bradt
While the nuclear shell model accurately describes the structure of nuclei near stability, the structure of
unstable, neutron-rich nuclei is still an area of active research. One region of interest is the set of nuclei
near N = 28. The shell model suggests that these nuclei should be approximately spherical due to the
shell gap predicted by their magic number of neutrons; however, experiments have shown that the nuclei
in this region rapidly become deformed as protons are removed from the spherical 48 Ca. This makes 46 Ar
a particularly interesting system as it lies in a transition region between 48 Ca and lighter isotones that are
known to be deformed.
An experiment was performed at the National Superconducting Cyclotron Laboratory (NSCL) to
measure resonant proton scattering on

46

Ar. The resonances observed in this reaction correspond to

unbound levels in the 47 K intermediate state nucleus which are isobaric analogues of states in the 47 Ar
nucleus. By measuring the spectroscopic factors of these states in 47 Ar, we gain information about the
single-particle structure of this system, which is directly related to the size of the N = 28 shell gap. Four
resonances were observed: one corresponding to the ground state in

47

Ar, one corresponding its first

excited 1/2− state, and two corresponding to 1/2+ states in either 47 Ar or the intermediate state nucleus.
However, only a limited amount of information about these states could be recovered due to the low
experimental statistics and limited angular resolution caused by pileup rejection and the inability to accurately reconstruct the beam particle track.
In addition to the nuclear physics motivations, this experiment served as the radioactive beam commissioning for the Active-Target Time Projection Chamber (AT-TPC). The AT-TPC is a new gas-filled
charged particle detector built at the NSCL to measure low-energy radioactive beams from the ReA3
facility. Since the gas inside the detector serves as both the tracking medium and the scattering target,
reactions are measured over a continuous range of energies with near-4π solid angle coverage. This experiment demonstrated that tracks recorded by the AT-TPC can be reconstructed to a good resolution,
and it established the feasibility of performing similar experiments with this detector in the future.

Copyright by
JOSHUA WILLIAM BRADT
2017

ACKNOWLEDGEMENTS

Although claiming that this work did not occur in a vacuum would be, at best, abuse of a metaphor, it is
certainly true that I could not have completed this thesis, let alone the years of work that came before
it, without the help of many others. First and foremost, I’d like to thank my advisor, Daniel Bazin, for his
endless patience and encouragement, and for giving me the freedom to experiment and solve problems
in my own way. Additionally, I am grateful to the entire AT-TPC group, which has made me feel welcome
from my first day; I truly could not have asked for a better group of people to work with. In particular,
I’d like to thank Wolfgang Mittig, Yassid Ayyad, and Michelle Kuchera for their help, feedback, and advice
throughout what sometimes seemed like a long and hopeless analysis process.
It would quite literally be impossible to finish graduate school without the help of my guidance committee, so I would like to thank Daniel Bazin, Morten Hjorth-Jensen, Wolfgang Mittig, Kirsten Tollefson,
and Remco Zegers for their advice and for taking the time out of their busy schedules to serve on my
committee. Additionally, I’d like to thank the staff of the NSCL for their role in this experiment, and the
MSU Institute for Cyber-Enabled Research for providing computing resources and technical support.
This whole experience would have been a lot more difficult and certainly a lot less fun without the
kind, supportive, and welcoming body of graduate students at MSU. I will truly miss all of you, but there
are a few people in particular that I’d like to mention. I’d like to thank Wei Jia Ong for conversations both
deep and superficial, for many cups of coffee, and for tolerating my puns. Also, thanks to Amy Lovell
for commiserating about Michigan drivers, providing baking advice, and playing no small part in me
passing my subject exams. I’d like to thank my officemates for tolerating the choice words I’ve addressed
to my computer many times in the last four years, and particularly I am grateful to Chris Izzo and Charlie
Loelius for many intriguing conversations about everything from pop culture to politics to geography.
Finally, I must thank the other graduate students from the AT-TPC group, including Lisa Carpenter and
Nate Watwood, for their company and for keeping me sane through long meetings and even longer shifts.
Last but certainly not least, I thank my parents for their endless love and support throughout my
life, for encouraging me to always do my best, and for helping me keep things in perspective when I
struggled. I know that wherever life takes me from here, I will always be able to count on your support,
and that means more than I could ever put into words.

iv

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

LIST OF ALGORITHMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix
CHAPTER 1
INTRODUCTION . . . . . . . . . . . . . . . . . . .
1.1 Nuclear physics background . . . . . . . . . . . . . . . . . .
1.1.1 Magic numbers and the nuclear shell model . . . .
1.1.2 Single-particle states . . . . . . . . . . . . . . . . . .
1.1.3 Isobaric analogue states . . . . . . . . . . . . . . . .
1.1.4 Elastic scattering in quantum mechanics . . . . . .
1.1.5 The optical potential . . . . . . . . . . . . . . . . . .
1.1.6 Resonances and the R matrix theory . . . . . . . .
1.1.7 Isobaric analogue resonances . . . . . . . . . . . .
1.1.8 Calculating the resonance energy . . . . . . . . . .
1.2 Experimental motivation . . . . . . . . . . . . . . . . . . . .
1.2.1 The disappearing N = 28 shell gap . . . . . . . . . .
1.2.2 Previous studies of 47 Ar . . . . . . . . . . . . . . . .
1.2.3 Resonant proton scattering in inverse kinematics .
1.2.4 Commissioning of the AT-TPC . . . . . . . . . . . .
1.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

CHAPTER 2
THE ACTIVE-TARGET TIME PROJECTION CHAMBER
2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Time projection chambers . . . . . . . . . . . . . . . . .
2.1.2 Active targets . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Structure of the AT-TPC . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Inner and outer gas volumes . . . . . . . . . . . . . . . .
2.2.2 Electric and magnetic fields . . . . . . . . . . . . . . . .
2.2.3 Sensor plane . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Tilting the detector . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Electronics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Data acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 Structure of the DAQ system . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

. 19
. 19
. 19
. 20
. 21
. 21
. 22
. 23
. 25
. 26
. 27
. 27

CHAPTER 3
ANALYZING AT-TPC DATA . . . . . . . . . . . .
3.1 Baseline correction . . . . . . . . . . . . . . . . . . . . .
3.2 Track reconstruction . . . . . . . . . . . . . . . . . . . .
3.3 Noise removal . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Hough transform for lines . . . . . . . . . . . . .
3.3.2 Hough transform for circles . . . . . . . . . . . .
3.3.3 Removing correlated noise . . . . . . . . . . . .
3.3.4 Removing uncorrelated noise . . . . . . . . . .
3.4 Modeling tracks and fits . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

. 31
. 31
. 32
. 35
. 35
. 38
. 39
. 40
. 43

v

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

1
1
1
2
4
5
6
7
9
11
11
11
13
16
17
18

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

CHAPTER 4
SIMULATIONS . . . . . . . . . . . . . . . . . . . . . . .
4.1 Overview of the simulation process . . . . . . . . . . . . . . . .
4.1.1 Parameter generation . . . . . . . . . . . . . . . . . . .
4.1.1.1
Uniform distribution . . . . . . . . . . . . . .
4.1.1.2
Simulated scattering distribution . . . . . .
4.1.2 Tracking and event generation . . . . . . . . . . . . . .
4.1.3 Simulating the trigger . . . . . . . . . . . . . . . . . . .
4.1.3.1
Generating channel-level trigger signals . .
4.1.3.2
Generating CoBo-level multiplicity signals .
4.1.3.3
Generating global triggers . . . . . . . . . . .
4.1.4 Cleaning and minimization . . . . . . . . . . . . . . . .
4.1.5 Post-processing and cuts . . . . . . . . . . . . . . . . .
4.2 Simulated acceptance . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Simulated resolution . . . . . . . . . . . . . . . . . . . . . . . .
4.4 R matrix simulation . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

. 58
. 58
. 58
. 59
. 61
. 64
. 65
. 66
. 67
. 68
. 68
. 68
. 69
. 86
. 89

CHAPTER 5
EXPERIMENT . . . . . . . . . . . . . . . . . . .
5.1 Beam production . . . . . . . . . . . . . . . . . . . . . .
5.2 Detection setup . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 AT-TPC setup . . . . . . . . . . . . . . . . . . . .
5.2.2 Ion chamber and secondary DAQ . . . . . . . .
5.2.3 Trigger . . . . . . . . . . . . . . . . . . . . . . . .
5.2.4 Excluding beam events . . . . . . . . . . . . . .
5.3 Analysis and results . . . . . . . . . . . . . . . . . . . . .
5.3.1 Correlating VME and GET events . . . . . . . .
5.3.2 Merging GET data files . . . . . . . . . . . . . . .
5.3.3 Beam particle identification . . . . . . . . . . .
5.3.4 Noise removal . . . . . . . . . . . . . . . . . . . .
5.3.5 Beam tracks . . . . . . . . . . . . . . . . . . . . .
5.3.6 Track fitting . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

3.5

3.6

3.7

3.4.1 Track representation and generation . . . . . . .
3.4.2 Simulated diffusion and track projection . . . . .
3.4.3 The objective function . . . . . . . . . . . . . . . .
3.4.4 Finding initial parameter values for optimization
The Monte Carlo optimizer . . . . . . . . . . . . . . . . .
3.5.1 Description of algorithm . . . . . . . . . . . . . .
3.5.2 Checking convergence . . . . . . . . . . . . . . . .
Alternative fitting algorithms . . . . . . . . . . . . . . . .
3.6.1 Gradient-based methods . . . . . . . . . . . . . .
3.6.2 Compass search . . . . . . . . . . . . . . . . . . .
3.6.3 Nelder-Mead simplex algorithm . . . . . . . . . .
3.6.4 Simulated annealing . . . . . . . . . . . . . . . . .
3.6.5 Neural networks . . . . . . . . . . . . . . . . . . .
3.6.6 Kalman filter . . . . . . . . . . . . . . . . . . . . .
Vertex energy reconstruction . . . . . . . . . . . . . . . .
3.7.1 Using the vertex position . . . . . . . . . . . . . .
3.7.2 Using kinematics . . . . . . . . . . . . . . . . . . .

vi

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

43
44
45
46
47
47
48
49
50
51
52
53
54
55
55
56
56

93
93
96
96
97
98
99
99
99
102
103
104
105
105

5.4

5.3.7 Proton parameter distributions . . . . .
5.3.8 Vertex energy reconstructions . . . . . .
5.3.9 Cross sections and excitation functions .
5.3.10 Alternative analysis method . . . . . . .
5.3.11 Results . . . . . . . . . . . . . . . . . . . .
5.3.12 Statistical testing . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . .
5.4.1 Comparison to previous measurements
5.4.2 Statistical significance . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

107
107
117
121
126
130
131
131
134

CHAPTER 6
CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
6.2 Future outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

vii

LIST OF TABLES

Table 1.1: First few levels in 47 Ar identified by Gaudefroy et al. [16, Table 1]. . . . . . . . . . . . . . .

14

Table 3.1: Monte Carlo fit parameters I , T , and R and initial parameter space widths σ used in
the 46 Ar(p, p ) measurement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

Table 4.1: Expected resonance properties calculated from the results of Gaudefroy et al. [18] using the procedure from Section 1.1.8. The resonance energy E R is given in MeV, and
the resonance width Γp , as calculated by DSIGMAIV, is given in keV. L is the orbital
angular momentum, J is the total angular momentum, Π is the parity, and S is the
spectroscopic factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

Table 4.2: Global optical model [29] parameters used in the simulation. Vr , R r , and a r are the
depth, radius, and diffuseness parameters for a real-valued Woods-Saxon potential,
while Vi , R i , and a i are the corresponding parameters for an imaginary Woods-Saxon
potential. The depth is given in MeV, and the radius and diffuseness are given in fm.
The Coulomb radius used was 1.27 fm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

Table 4.3: Trigger conditions used in the experiment. As noted in the text, the unusually large
number of excluded pads in set 2 was due to a user error. . . . . . . . . . . . . . . . . . .

69

Table 4.4: Results of fitting the distributions of parameter and energy deviations from Figs. 4.19
and 4.21. The Laplace and Gaussian distributions are defined in Eqs. (4.22) and (4.23),
respectively. The units listed apply to both the centroid and width parameters. The
full width at half-maximum (FWHM) was defined as 2σ 2 log(2) for the Gaussian
distributions and 2b log(2) for the Laplace distributions. . . . . . . . . . . . . . . . . . . .

87

Table 5.1: Parameter values used while analyzing the 46 Ar(p, p ) data. The vectorial values are
given in the beam coordinate system of Fig. 3.2. . . . . . . . . . . . . . . . . . . . . . . . . 106
Table 5.2: Properties of the resonances shown in Fig. 5.26. Energies and widths are given in
CM
keV/u. E Ar
is the energy of the corresponding level in the analogue nucleus calculated assuming ∆EC − S n = 2807(30) keV, from the literature values of those parameters; E x is the excitation energy calculated from the resonance assuming ∆EC − S n =
2680(20) keV, the energy of the ground state resonance; S refers to the spectroscopic
factor; Γ is the total resonance width; Γp is the proton single-particle width; and Γ↓
is the spreading width arising from mixing with the T< states. For quantities with two
uncertainties, the first value is the systematic uncertainty and the second is statistical. . 127
Table 5.3: Renormalization factors applied to spectroscopic factors from DSIGMAIV. . . . . . . . . 129
Table 5.4: Calculated values of the resonance integral from Eq. (5.17) and the relative uncertainty δI /I it contributes to the spectroscopic factor for a variation of 10 keV in Γ↓ . . . . 129
Table 5.5: Straggling at the resonance locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
viii

Table 5.6: Statistical significance of the proposed resonances. The first three lines show the results of the F -test applied to the three individual resonances. The last line shows the
result of the test applied to the data as a whole. . . . . . . . . . . . . . . . . . . . . . . . . 131
Table 5.7: Literature values of the measured excitation energies and spectroscopic factors. The
SDPF, SDPF-U, and SDPF-MU values are from shell model calculations performed by
the authors of the referenced papers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

ix

LIST OF FIGURES

Figure 1.1: Chart of the nuclides. Neutron number increases along the horizontal axis, and proton number increases along the vertical axis. The color corresponds to the log of the
half-life, with lighter colors indicating a longer half-life. Stable nuclei are colored black.

2

Figure 1.2: The first few levels of the nuclear shell model (level spacings are not drawn to scale).
The number in parentheses next to each level indicates the number of nucleons that
fit in that level. The shell gaps are shown with dotted lines. The boxed number to the
left of each shell gap line is the number of nucleons that fit below that shell closure.
These are the first few magic numbers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

Figure 1.3: Three isobaric nuclei with isospin projections T z = −1, 0, and +1 from left to right.
Sets of isobaric analogue states are indicated with dashed lines. The ground states of
the 6 He and 6 Be spectra are printed higher than that of 6 Li to make the correspondence more apparent. Energies are given in keV. (Data taken from [38].) . . . . . . . . .

5

Figure 1.4: A high-resolution measurement of the excitation function from 40 Ar(p, p ). The bottom plot shows a theoretical fit to the data near 2.45 MeV. Spin-parity assignments
were made based on the shape of the resonances, and it is apparent that the 1/2+ and
3/2− resonances are spit into several fine-structure components. (Figure reproduced
from [22, Fig. 19.13], which was adapted from [7].) . . . . . . . . . . . . . . . . . . . . . .

10

Figure 1.5: Calculating the resonance energy of the isobaric analogue state labeled “IAS”. In this
figure, S n is the neutron separation energy, δ is the difference between the neutron
and proton masses, ∆E c is the Coulomb energy of the last proton, and E r is the resonance energy. The ground state of 47 K is located far below the bottom of the figure
and is not shown. The energy separations are not drawn to scale. . . . . . . . . . . . . .

10

Figure 1.6: Trends in E (2+
1 ) and B (E 2) near the N = 20 and N = 28 magic numbers. (Reproduced
from [49].) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

Figure 1.7: Energy spectrum of d 46 Ar, 47 Ar p observed by Gaudefroy et al. [16]. Subplot (a) shows
the spectrum with background curves from carbon-induced reactions only (labeled
“C”) and from carbon-induced reactions and deuteron breakup (labeled “C + d”).
Subplot (b) shows the background-subtracted spectrum fit with 9 Gaussian curves.
Subplot (c) shows the spectrum of events where protons were in coincidence with
the 47 Ar nucleus. (Reproduced from [16, Fig. 1].) . . . . . . . . . . . . . . . . . . . . . . .

14

Figure 1.8: A comparison of neutron single-particle energies in 49 Ca and 47 Ar from Gaudefroy
et al. [16]. The splitting of the 1p orbitals is reduced in 47 Ar compared to 49 Ca. (Reproduced from [16, Fig. 3].) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

Figure 1.9: Energy levels identified by Bhattacharyya et al. [6]. (Adapted from [6, Fig. 3].) . . . . . .

16

x

Figure 1.10: γ-ray spectrum measured by Gade et al. [15]. The results from each reaction are
shown separately, and levels indicated in gray were not observed using a given reaction. (Reproduced from [15, Fig. 3].) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

Figure 2.1: Principle of operation of a TPC. A uniform electric field is created between the anode
and the cathode. This field transports the ionization electrons produced by a charged
particle towards the anode, where they are amplified and collected. . . . . . . . . . . .

19

Figure 2.2: A schematic view of the AT-TPC. The outer shielding volume was made transparent in this image to make the details of the inner volume more visible. Beam enters
the detector through the beam duct at the right-hand side of the image and moves
toward the sensor plane on the left. Some components of the GET electronics are
shown mounted on the downstream end of the detector (see Section 2.4). . . . . . . . .

21

Figure 2.3: A photograph of the inner volume wall with the rings of the inner field cage installed.
The outer field cage rings have a larger diameter, but are otherwise identical and are
installed on the outer surface of the cylinder. . . . . . . . . . . . . . . . . . . . . . . . . .

22

Figure 2.4: Layout of the pad plane. The inset shows a closer view of one corner of the hexagonal inner region. This region of half-scale pads provides finer resolution near the
reaction vertex, which will generally occur near the central axis of the detector. . . . . .

24

Figure 2.5: Principle of operation of a Micromegas. The electric field magnitudes shown are
nominal: the field in the region above the mesh may vary from roughly 103 V/m to
105 V/m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

Figure 2.6: Illustration of pileup separation due to the tilt angle. The leftmost track is elastic
α-α scattering, and the rightmost track is another α particle that entered the detector during the event. Without the tilt angle, the two beam tracks would be collinear,
making them more difficult to separate. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

Figure 2.7: A schematic view of the GET electronics system. Only one AsAd is shown for clarity,
but 40 AsAds are used to instrument the AT-TPC. . . . . . . . . . . . . . . . . . . . . . . .

26

Figure 2.8: A photograph of the fully instrumented AT-TPC mounted inside its solenoid magnet. The 40 AsAds are mounted on the downstream end flange surrounded by copper
shielding. Individual cables connect each AsAd to its controlling CoBo, and a separate set of cables connects each to its power supply. . . . . . . . . . . . . . . . . . . . . .

28

Figure 2.9: Structure of the DAQ system. Only two worker nodes are shown, but there are ten in
the actual system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

Figure 3.1: Subtracting the baseline from a single signal. The top plot shows the raw signal and
the calculated baseline. The bottom plot shows the difference between the two. . . . .

32

xi

Figure 3.2: Coordinate systems used in the analysis. After calibration, the recorded tracks are
in the detector coordinate system (a). This system has its origin in the center of the
sensor plane at z = 0 with the z axis pointing upstream toward the beam entrance
window. The beam coordinate system (b) is related to the detector coordinate system by a rotation through angle τ about the beam entrance window, where τ is the
detector tilt angle. In this system, the beam moves along the −w axis. Note that both
of these coordinate systems are left-handed. . . . . . . . . . . . . . . . . . . . . . . . . .

33

Figure 3.3: Event 99-30 from the 46 Ar data set, a particularly noisy event, before cleaning. There
are some correlated noise structures present at the bottom of the uv projection and
the left of the w v projection, and there is uncorrelated noise throughout the event. . .

36

Figure 3.4: Parameterizing a straight line in two dimensions using a normal line drawn to the origin. 37
Figure 3.5: Transforming points into the Hough space. The left-hand plot shows three points
A, B , and C in the normal coordinate space with lines drawn between them. These
points map to sinusoidal curves in the Hough space, as shown in the right-hand plot.
Lines in the normal system become points in the Hough space, so a line that passes
through two points in the normal space is the intersection of two curves in the Hough
space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

Figure 3.6: Finding the center of a circle using the Hough transform. . . . . . . . . . . . . . . . . . .

38

Figure 3.7: Cleaning an event with the Hough transforms. Plot (a) shows the data from Fig. 3.3
in the ρφ vs w space with the found lines. The colored points are classified as data,
while the empty points are noise. Plot (b) shows the Hough space calculated for this
data. The maximum angular slice is shown between the two white lines, and the 1D
projection of this slice is shown on the right with the peak locations marked. . . . . . .

41

Figure 3.8: Final results of the noise removal process on the event shown in Fig. 3.3. The blue
points were kept, and the semi-transparent red points were removed as noise. Some
noise is still present after cleaning, but much of it was removed. . . . . . . . . . . . . . .

42

Figure 3.9: Template used to simulate diffusion. The charge q from the initial point is distributed
among the diffused points. The central point gets 0.4q units of charge, and the outer
points each get 0.6q/8 units of charge. This is based on a two-dimensional normal
distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

Figure 3.10: Samples drawn by the Monte Carlo algorithm while fitting one event. The horizontal
axes represent the sample number, while the vertical axes represent the value of each
of the 6 variables. x, y, and z are given in meters, E is in MeV, and θ and φ are in
radians. Each 500-sample iteration of the code is visible as a large-scale block in this
plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

xii

Figure 3.11: Results of the Monte Carlo convergence study. Subplots represent different numbers
of points T per iteration. Within each plot, rows represent the reduction factor R, and
columns represent the number of iterations I . The number in each cell represents
the percentage of fits in which the optimizer converged to the correct minimum,
subject to a tolerance of 10 mm in each dimension for the vertex position, 100 keV for
the initial proton energy, 5° for the azimuthal angle, and 2° for the scattering angle.
The fit was run 1000 times on the same simulated track for each set of parameters. . .

49

Figure 3.12: Simulated two-dimensional error surface and a one-dimensional slice through that
surface. These were generated from a simulated 46 Ar(p, p ) scattering event by varying the scattered proton’s initial energy and scattering angle, calculating the objective function at each point. 400 points were simulated in each dimension. On the
surface plot, the horizontal axis represents energy in MeV, the vertical axis represents the energy component of the objective function, and the axis going into the
page represents the scattering angle in radians. . . . . . . . . . . . . . . . . . . . . . . . .

50

Figure 3.13: An illustration of the points tested in the Nelder-Mead simplex method. The vertex v 2
is reflected across the point c using four different values of the multiplicative factor
α, leading to the four possible new vertices, shown in purple, yellow, blue, and green.
(Reproduced from [28, Fig. 4.2(b)].) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

Figure 3.14: Steps taken by a vertex during simplex collapse, a failure mode of the Nelder-Mead
algorithm. (Reproduced from [35, Fig. 2.1].) . . . . . . . . . . . . . . . . . . . . . . . . . .

53

Figure 4.1: Illustration of coordinate calculation for uniform distributions. The meanings of the
variables are given in Section 4.1.1.1. Note that this is a projection into the z y plane,
and the true three-dimensional orientation also depends on azimuthal angles φ0
and φB and vertex location x 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

Figure 4.2: Results of DSIGMAIV calculation in the region of interest. . . . . . . . . . . . . . . . . . .

62

Figure 4.3: Schematic of a channel in the GET electronics, showing the trigger threshold circuit.
(Reproduced from [3].) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

Figure 4.4: Trigger exclusion and low-gain regions used in the simulation and in the experiment.
Only pads in the light and dark gray regions contribute to the trigger. . . . . . . . . . . .

70

Figure 4.5: Distributions of the χ2 values for the simulated events in the first three sets of runs.
Events that fell outside of the shaded region between the two vertical black lines
on the histogram were thrown away. The events that were kept had χ2pos < 1.0 and
χ2en < 1.0. This plot does not include events for which the fitter returned an error. . . .
Figure 4.6: Distributions of χ2 values for the final two sets of runs. See the caption of Fig. 4.5 for
details. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xiii

71

72

Figure 4.7: Simulated acceptance of the trigger and analysis as a function of the track parameters for trigger parameter set 0. The trigger and analysis cuts were applied independently to the simulated events, leading to six classes of events. Events are grouped
into “trigger” and “no trigger” groups depending on whether they would have triggered the electronics in a real experiment. These groups are then subdivided into
“good fit” and “bad fit” groups depending on whether they pass the χ2 cuts. Events
in the “fit failed” groups caused an error in the Monte Carlo algorithm. Only events in
the “trigger, good fit” group would have been recorded and correctly reconstructed
in the experiment. The six groups are stacked in these histograms, so the contour of
the top of each plot reflects the shape of the total distribution. . . . . . . . . . . . . . . .

73

Figure 4.8: Simulated acceptance of the trigger and analysis for trigger parameter set 2. The large
trigger exclusion region used in these runs prevented events outside of a very small
angular and energy domain from being accepted. . . . . . . . . . . . . . . . . . . . . . .

74

Figure 4.9: Simulated acceptance of the trigger and analysis for trigger parameter set 3. . . . . . .

75

Figure 4.10: Simulated acceptance of the trigger and analysis for trigger parameter set 4. . . . . . .

76

Figure 4.11: Simulated acceptance of the trigger and analysis for trigger parameter set 5. The acceptance was reduced for these runs due to the exclusion of CoBos 4 and 9 from the
trigger. This is most visible in the distribution with respect to the azimuthal angle. . . .

77

Figure 4.12: Efficiencies of the AT-TPC trigger, the analysis χ2 cuts, and the total of the two for the
first three sets of parameters, plotted as functions of the 46 Ar energy at the reaction
vertex and the scattering angle in the center of mass frame. The bins used were 5°
wide in the angular dimension and 50 keV wide in the energy dimension. . . . . . . . .

79

Figure 4.13: Efficiencies of the AT-TPC trigger, the analysis χ2 cuts, and the total of the two for the
last two sets of parameters. See the caption of Fig. 4.12 for details. . . . . . . . . . . . .

80

Figure 4.14: Efficiencies for parameter set 0 as a function of 46 Ar vertex energy at several scattering angles. The error bars shown are ±σ ,i as defined in Eq. (4.21). The bins used in
this plot are the same as those used in the 2D histograms in Figs. 4.12 and 4.13. . . . .

81

Figure 4.15: Efficiencies for parameter set 2 as a function of 46 Ar vertex energy at several scattering angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

Figure 4.16: Efficiencies for parameter set 3 as a function of 46 Ar vertex energy at several scattering angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

Figure 4.17: Efficiencies for parameter set 4 as a function of 46 Ar vertex energy at several scattering angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

Figure 4.18: Efficiencies for parameter set 5 as a function of 46 Ar vertex energy at several scattering angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85

xiv

Figure 4.19: Simulated resolution of the AT-TPC. These six plots show the error in the reconstruction of each of the six track parameters, defined as the reconstructed value minus
the correct value. The results of the fits to these distributions are given in Table 4.4. . .

87

Figure 4.20: Illustration of relationship between uncertainty in vertex location and error in azimuthal angle. The reconstructed vertex is marked with an orange square, while the
actual vertex location is marked with a blue circle. The red curve represents the proton track. The tangent lines at these two locations have very different slopes. . . . . . .

88

Figure 4.21: Simulated resolution of the two energy reconstruction techniques from Section 3.7.
The widths of the distributions are given in Table 4.4. ∆E kin and ∆E pos are defined
as the difference between the energies as calculated from the reconstructed parameters and the known correct parameters. For the rightmost plot of E kin − E pos , both
values were found from the reconstructed parameters. The difference between the
two energies for the known correct parameters is not shown since those two values
are identical. Note that the horizontal scale is different on the center plot. . . . . . . . .

89

Figure 4.22: Distributions of the objective function components for simulated events drawn from
the R matrix calculation results. Each distribution was cut at 1.0, discarding all events
with a larger value of χ2 . This is the same cut that was used for the flat distribution
(see Figs. 4.5 and 4.6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

Figure 4.23: Two-dimensional histogram of results from the R matrix–based simulation. The two
expected resonance locations are marked on the vertical axis at 2.696 MeV/u and
3.875 MeV/u. The two individual bright bins at 80° and 85° are artifacts from dividing
by the very small number of counts in those bins. . . . . . . . . . . . . . . . . . . . . . .

90

Figure 4.24: Excitation functions derived from the R matrix simulation. The expected locations
of the two resonances are shown with vertical orange lines. . . . . . . . . . . . . . . . .

92

Figure 5.1: A diagram of the NSCL beamline, with some regions of the facility labeled. The beam
starts at the ion sources (label 1) and ends at the AT-TPC (label 14). (Reproduced
from http://nscl.msu.edu/public/virtual-tour.html.) . . . . . . . . . . . . . . .

94

Figure 5.2: The A1900 fragment separator and the coupled cyclotrons. (Reproduced from [36].) . .

94

Figure 5.3: Timing of events recorded by the GET electronics. Particles arrived in bunches at a
rate of approximately 2 Hz. Note that this only shows events that triggered the electronics; in reality, many more nuclei arrived in each bunch. . . . . . . . . . . . . . . . .

95

Figure 5.4: A signal from the ion chamber (see Section 5.2.2) from an event where four particles
entered the detector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

xv

Figure 5.5: Plots of various signals in a spark-like event. The top left plot shows the individual
signals in the GET electronics, and the top right plot shows the sum of these signals.
The bottom left plot contains the signal from the ion chamber, and the bottom right
plot shows the signal from the Micromegas mesh. The gray regions indicate domains
in which there is no data. The time buckets in the two DAQ systems are of the same
width, but are offset by 100 due to the trigger delay. . . . . . . . . . . . . . . . . . . . . . 100
Figure 5.6: Plots of (x, y, z) positions corresponding to peak locations in a spark-like event. The
majority of the channels in the inner region of the pad plane fired at approximately
the same time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Figure 5.7: A raw ion chamber signal produced by an 46 Ar nucleus. The passage of the nucleus
through the chamber is indicated by the large negative peak around t = 160. The
ADC used to record this signal was 12 bits, but it also had an overflow bit that made
it effectively record 13 bits per sample for a maximum digitized value of 8192. . . . . . 103
Figure 5.8: Distribution of peak heights in the ion chamber. The center of the 46 Ar peak is indicated by the dashed vertical line, and the boundaries of the cut applied around this
peak are shown with solid lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
Figure 5.9: A recorded beam track from the 46 Ar experiment. The top panel shows the traces
from the low-gain region of the sensor plane, which is where the beam track is projected onto the pads. The signals have an unusual shape which is much different
from the signals generated by the proton tracks. The inset illustrates the failure of
the baseline restoration for one track. In the lower panel, the red squares correspond
to the reconstructed beam track, which is quite sparse. . . . . . . . . . . . . . . . . . . . 109
Figure 5.10: Distributions of the objective function components from fitting experimental data.
The position component was cut at 14, while the energy component was cut at 10.
Events with larger values for either of these components were discarded. . . . . . . . . 110
Figure 5.11: Distribution of reconstructed proton vertex energies from the experiment. Note that
the vertical axis uses a logarithmic scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
Figure 5.12: Distribution of reconstructed scattering angles in the center-of-mass frame. . . . . . . 111
Figure 5.13: Distribution of vertex positions in the detector. Plot (a) shows the overall distribution
of vertex positions in the plane transverse to the beam. Plot (b) shows that the radius
of this distribution increases for smaller z, or for vertex positions further from the
beam entrance window. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
Figure 5.14:

46

Ar vertex energies as reconstructed from the vertex position (top) and from kinematics (bottom). The reconstruction methods are detailed in Section 3.7. . . . . . . . . 113

Figure 5.15: A comparison of the two vertex energy reconstructions. The line y = x is plotted for
reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

xvi

Figure 5.16: Positional vertex energy recalibration. The 2D histogram shows values of E p / cos2 (θlab ),
and the two curves show values of 4E Ar . The original calibration is shown with a
dashed red curve, and the corrected calibration is shown with a solid cyan curve. The
error band on the corrected calibration shows the interval 1.065 ≤ α ≤ 1.115, where
the best value is in the center at α = 1.09. Events with θlab > 65° were eliminated from
this plot to reduce the spread at higher energies caused by cos2 (θlab ) approaching 0. . 116
Figure 5.17: Derivatives of kinematic lines at a few proton energies. The vertex energy changes
very rapidly as a function of scattering angle as the scattering angle approaches 90°
in the laboratory frame. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
Figure 5.18: Results of a fit to the edge of the positional vertex energy distribution. The fit was performed with the function defined in Eq. (5.4), yielding parameters µ = 4.11 MeV/u,
σ = 19.7 keV/u, A = 694, and b = 1.98. The bins in the energy dimension were 10 keV/u
wide, and the error bars shown are ± N . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Figure 5.19: Two-dimensional histogram showing measured counts binned as a function of 46 Ar
vertex energy in the laboratory frame and scattering angle in the center-of-mass frame. 118
Figure 5.20: Unnormalized excitation functions as a function of 46 Ar vertex energy in the laboratory frame, shown at a variety of center-of-mass scattering angles. The error bars
shown are ± N , where N is the number of counts in a given bin. . . . . . . . . . . . . . 119
Figure 5.21: Simulated total efficiencies for parameter set 3 (as defined in Table 4.3) with kernel
ridge regression fit. Similar fits were performed for the efficiencies generated for the
other parameter sets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
Figure 5.22: Excitation functions in the region of interest for 46 Ar(p, p ) calculated from experimental results. The results of an R matrix calculation performed with DSIGMAIV is
also shown for comparison. The parameters used in this calculation are given in Section 5.3.10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
Figure 5.23: Experimental counts summed over all scattering angles and fit with a quadratic function. The fit function was y(x) = −263(12)x 2 +1687(74)x −1207(108), where the numbers in parentheses are the standard errors of the fit coefficients. A 1σ error band is
shown around the fit using shading. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
Figure 5.24: Weighted total excitation function defined in Eq. (5.11) calculated for both the full R
matrix simulation (labeled “Full calculation”) and a calculation with no resonances
(labeled “Baseline”). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
Figure 5.25: Effect of convolving the R matrix calculation with a Gaussian. The resonances in the
convolved version are wider and have significantly smaller amplitudes. . . . . . . . . . 126
Figure 5.26: Comparison of data to R matrix calculation. . . . . . . . . . . . . . . . . . . . . . . . . . 127

xvii

Figure 5.27: Results of energy straggling calculation with TRIM [57]. The orange line is a linear
fit to the data with slope −3.12(35) keV/MeV and intercept 16.2(10) keV. The shaded
orange band is a 1σ error band. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
Figure 5.28: Comparison between experimental results and literature values. Filled shapes indicate experimentally measured values, while empty shapes indicate shell model calculations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

xviii

LIST OF ALGORITHMS

Algorithm 1:

Hough space binning method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

Algorithm 2:

Counting nearest neighbors inside a radius . . . . . . . . . . . . . . . . . . . . . . . .

42

Algorithm 3:

Basic compass search algorithm [28] . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

Algorithm 4:

Basic simulated annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

xix

CHAPTER 1
INTRODUCTION

One of the primary goals of nuclear physics research is to understand the structure of the nucleus. While
the introduction of the nuclear shell model helped explain the structure of stable nuclei, the exotic,
neutron-rich nuclei produced by modern radioactive beam facilities continually present new challenges
to our understanding of nuclear structure far from stability. In these neutron-rich nuclei, the predictions
made by the shell model based on stable nuclei begin to break down, leading to changes in single-particle
levels and shell gaps.
This thesis describes an experiment that was performed to study the structure of

47

Ar, an isotope

with 18 protons and 29 neutrons that lies in a region of particular interest for studies of shell evolution.
This introductory chapter begins with a brief review of some relevant concepts from nuclear physics.
The motivations behind the experiment are then described in detail.

1.1 Nuclear physics background
Individual nuclei are characterized by the number of neutrons and protons they contain. A group of
nuclei with the same number of protons Z are called isotopes, while a group with the same number of
neutrons N are called isotones. Nuclei with the same mass A = Z + N are called isobars. The full set of
known nuclei is often displayed in a chart like the one in Fig. 1.1. There, vertical lines group isotones,
horizontal lines group isotopes, and lines parallel to the line Z = −N group isobars.
1.1.1 Magic numbers and the nuclear shell model
Although several models have been proposed to describe the structure of the nucleus, the best model
currently known is the nuclear shell model. The shell model was first developed in 1949 [9] by physicists
including Maria Goeppert Mayer, J. Hans D. Jensen, and Eugene Wigner—who later shared the 1963 Nobel Prize in Physics for their work—to explain the existence of the so-called magic numbers of protons
and neutrons shared by many of the most stable nuclei. The key feature of this model was the proposal
of a shell structure for nucleons (protons and neutrons) analogous to that of electrons in atomic orbitals.

1

120

Number of Protons

100

80

60

40

20

0
0

50

100

150

Number of Neutrons

Figure 1.1: Chart of the nuclides. Neutron number increases along the horizontal axis, and proton number increases along the vertical axis. The color corresponds to the log of the half-life, with lighter colors
indicating a longer half-life. Stable nuclei are colored black.
A few of the lowest-energy shells are shown in Fig. 1.2. The left-hand column shows the shells predicted by assuming that nucleons are bound by a harmonic oscillator potential. However, to reproduce
many of the experimentally observed magic numbers (particularly those at higher N or Z ), the model
must be extended to include a spin-orbit interaction [49]. This produces the set of levels shown on the
right-hand side of the figure. For example, the neutron shell gap at N = 28 arises due to spin-orbit splitting of the 1 f oscillator shell that pushes the 1 f 7/2 level to a somewhat lower energy than the 1 f 5/2 level
and the two levels derived from the 1p shell, which are all grouped together at a higher energy.

1.1.2 Single-particle states
The simple picture of the shell model presented above is derived for a single particle bound in a simple
potential well, like that of a harmonic oscillator or Woods–Saxon function. However, the energy levels of
actual nucleons are much more complex and depend on interactions with all of the other nucleons in the
system. Therefore, the energy levels measured experimentally are not equivalent to the single-particle
energy levels predicted by theory. Instead, the single-particle energies can be defined as a weighted sum

2

1g
50
1g 9/2 (10)
2p 1/2 (2)
2p

28

1 f 5/2 (6)
2p 3/2 (4)

1f

1 f 7/2 (8)
20
2s
1d

1d 3/2 (4)
2s 1/2 (2)
1d 5/2 (6)

1p

1p 1/2 (2)
1p 3/2 (4)

1s

1s 1/2 (2)

8

2

Figure 1.2: The first few levels of the nuclear shell model (level spacings are not drawn to scale). The
number in parentheses next to each level indicates the number of nucleons that fit in that level. The
shell gaps are shown with dotted lines. The boxed number to the left of each shell gap line is the number
of nucleons that fit below that shell closure. These are the first few magic numbers.
of the actual energy levels: [2]
α

In this sum, the single-particle energy

α

=

S(α, i )E i .

(1.1)

i

of particle α is related to the experimental energy levels E i

having the same spin and parity as the single-particle state. The weighting factors S(α, i ) are the spectroscopic factors, which can be thought of as the probability of finding the particle α in the given singleparticle state.
A more thorough theoretical treatment of the notion introduced in Eq. (1.1) is given by Baranger [2].
They define a single-nucleon Hamiltonian W by
〈α|W |β〉 =

n

〈ψ0 |a α |χn 〉 (E n − E 0 ) 〈χn |a β† |ψ0 〉 +

N

〈ψ0 |a β† |χN 〉 (E 0 − E N ) 〈χN |a α |ψ0 〉 .

(1.2)

The matrix elements in this sum connect the ground state |ψ0 〉 (with energy E 0 ) of a nucleus with mass
number A to the set of states |χn 〉 (with energy E n ) in the adjacent A + 1 “particle” nucleus and the set
3

of states |χN 〉 (with energy E N ) in the A − 1 “hole” nucleus. |β〉 and |α〉 are two single-particle states with
corresponding creation and annihilation operators a β† and a α . Baranger [2] points out that if these states
are assumed to be eigenstates of W , then W must be diagonal, so Eq. (1.2) reduces to

α

≡ 〈α|W |α〉 =

n

| 〈ψ0 |a α |χn 〉 |2 (E n − E 0 ) +

N

| 〈χN |a α |ψ0 〉 |2 (E 0 − E N ).

(1.3)

They then identify | 〈ψ0 |a α |χn 〉 |2 with the spectroscopic factor and (E n − E 0 ) with the energy E i from
above, assuming the reference state |ψ0 〉 is the ground state.
1.1.3 Isobaric analogue states
Many measurements of nuclear spectra have shown that families of isobaric nuclei, or nuclei with the
same mass number A but different charges Z , tend to share sets of common energy levels. This occurs
because the neutron and proton have very similar masses and behave approximately identically under
the strong nuclear force, leaving charge as the main difference between them. If Coulomb repulsion is
neglected, the proton and neutron can then be considered as two states of the same particle.
Formally, this symmetry is described using the isospin t, a vector quantity that behaves identically to
ordinary spin vectors. Protons are defined1 as having isospin projection t z = −1/2, while neutrons have
t z = +1/2. The total isospin T of a nucleus is the vector sum of the isospins of its constituent nucleons,
and the total isospin projection, which distinguishes individual members of the isobaric multiplet, can
be altered by defining a set of isospin raising and lowering operators much like those used for spin and
orbital angular momentum.
Much like how a magnetic field splits degenerate electron orbitals into distinguishable levels based
on the projection of their total angular momentum, the presence of the Coulomb force inside the nucleus
creates multiplets of states with the same total isospin but different isospin projections. The rest of the
quantum numbers of these states are the same, so they have very similar wave functions.
Figure 1.3 shows an example of isobaric analogue states in three nuclei with A = 6. Two isospin
triplets (T = 1) are visible, connecting sets of 0+ and 2+ states. The remaining, unique states in 6 Li are
isospin singlets (T = 0).
1 Although this sign convention is commonly used in nuclear physics, other fields use the opposite convention.

4

2+

1797

0+

0
6

2+

5366

2+

4312

0+

3563

(2)+

1670

0+

0
6

He
3+

Be

2186

1+

0
6

Li

Figure 1.3: Three isobaric nuclei with isospin projections T z = −1, 0, and +1 from left to right. Sets of
isobaric analogue states are indicated with dashed lines. The ground states of the 6 He and 6 Be spectra
are printed higher than that of 6 Li to make the correspondence more apparent. Energies are given in keV.
(Data taken from [38].)
1.1.4 Elastic scattering in quantum mechanics
The quantum mechanical description of the elastic scattering of nuclei off a potential V (r) depends on
the solution of the time-independent Schrödinger equation [22]
2m

∇2 ψ +

2

[E − V (r)] ψ = 0.

(1.4)

In spherical coordinates, the solution to this equation must approach
r →∞

ψ(r) −−−−→ e i kz +

e i kr
f (θ)
r

(1.5)

at large distances from the scattering center at the origin [22]. The first component of this expression
represents an unscattered plane wave with wavenumber k moving along the same direction z as the
incoming particle, and the second component represents a spherical wave propagating outward from
the scattering center as a function of radius r . Assuming azimuthal symmetry, the factor f (θ) defines
what proportion of the incoming wave was scattered as a function of the scattering angle θ, so the elastic
scattering cross section can be written as [22]
dσ
2
= f (θ) .
dΩ

(1.6)

Equation (1.4) is commonly solved by expanding ψ over partial waves with discrete values of angular

5

momentum L:
ψ(r, θ) =

∞

u L (r )
P L (cos θ).
r
L=0

(1.7)

Here, u L (r ) is the radial component of the wave function, and the P L (cos θ) are the Legendre polynomials. Applying boundary conditions leads to the following solution for f (θ): [22]
f (θ) =

1
2i k

∞
L=0

(2L + 1) e 2i δL − 1 P L (cos θ).

(1.8)

Here, k is the wavenumber, and δL is called the phase shift, which contains all of the physical information
about the scattering [22]. The phase shift is often used to define the scattering or S matrix element
S L = e 2i δL .

(1.9)

1.1.5 The optical potential
One commonly used model for the nuclear potential is the optical model. This model represents the nuclear potential as a complex quantity with a real component that produces scattering and an imaginary
component that absorbs incoming flux [22]. The absorptive component accounts for non-elastic processes in the interaction [22]. This is known as an “optical” model since it was created by analogy to the
scattering of light in a semi-opaque medium, where the proportion of light absorbed by the medium can
be quantified using the imaginary part of a complex index of refraction [22].
An optical model is typically formulated using the Woods-Saxon potential
V (r ) = −

V0
1 + exp

r −R
a

,

(1.10)

where V0 is the depth of the potential well, R is its radius, and a is its diffuseness, a measure of how quickly
the function rises [51]. Both the real and imaginary parts of the optical potential have the form shown in
Eq. (1.10), but they may have different parameters.
The nuclear potential is parameterized in terms of several Woods-Saxon potentials. There is generally
a real and imaginary volume potential term, and an imaginary surface term is often included as well [51].
Finally, a spin-orbit term may be added to account for the coupling between the spin of the nucleus and
its angular momentum [51].
The potential described above only contains a contribution from the strong nuclear force, which acts
at short ranges. Scattered charged particles, such as protons, will also experience a long-range Coulomb
6

force due to the nuclear charge. This contributes the following term to the total potential: [51]

2

VCoul. (r ) = Zproj. Z e ×






3
2

2

− 2Rr2



1

Coul.

1
R Coul.

r ≤ R Coul.

(1.11)

r ≥ R Coul.

r

This potential term depends on the charge Zproj. e of the projectile, the charge Z e of the nucleus, and the
Coulomb radius R Coul. , which is typically proportional to A 1/3 in a nucleus with mass number A [51].
In total, this leaves us with a model of the scattering potential that depends on 12 or more free parameters which are found by fitting the model to experimental data. Fitting the model to data from a single
nucleus at a single energy produces a local optical potential, and these local parameterizations will generally describe the data best [51]. However, the model can also be fit to data from a range of nearby nuclei
over a range of energies to produce a global optical potential with parameters that vary slowly as functions of A and energy. These global potentials tend to fit scattering data worse than a local potential for
the nucleus being studied, but they have the benefit of allowing interpolation or extrapolation to nuclei
for which a reliable set of data does not yet exist [51].

1.1.6 Resonances and the R matrix theory
During a low-energy nuclear reaction, it is possible for the projectile and target nuclei to temporarily
form a short-lived, intermediate-state compound nucleus whose decay products correspond to the exit
channel of the reaction. This compound nucleus may be formed in any of a large number of unbound
energy levels which correspond to resonances in the measured excitation function, or isolated regions
of the excitation function where the cross section departs significantly from the baseline value. As the
density of states increases with increasing excitation energies, these resonances eventually merge into a
continuous distribution at high energies [22].
Resonances typically occur in one reaction channel and one partial wave, giving them well-defined
values of orbital angular momentum L and total angular momentum J . For a spin-0 particle, the cross
section near an isolated resonance of this type at energy E 0 with width Γ is described with the Breit–
Wigner formula [22]
(2L + 1)2
Γ2
dσ
(E , θ) =
P 2 (cos θ),
dΩ
4k 2 (E − E 0 )2 + 14 Γ2 L
as a function of energy E and scattering angle θ, with wavenumber k.

7

(1.12)

Resonances in measured elastic scattering cross sections can be fit using the phenomenological R
matrix theory. In the formulation used by Descouvemont and Baye [11], the R matrix for an isolated
resonance in a single partial wave with orbital angular momentum
R (E ) =

γ2r
Er − E

is given by

.

(1.13)

In this equation, E r and γr are the formal resonance energy and reduced width. These can be calculated
from the observed resonance energy E R and width ΓR as follows. The formal reduced width γr can be
calculated from the observed reduced width γR using the equation
γ2r =

γ2R
1 − γ2R S (E R )

.

(1.14)

The observed reduced width is found from the observed resonance width ΓR :
γ2R =

ΓR
.
2P (E R )

(1.15)

The formal resonance energy E r can be calculated in a similar way:
E r = E R + γ2r S (E R ).

(1.16)

In Eqs. (1.14) to (1.16) above, S is called the shift factor, and S is its derivative. P is the penetration factor
and should not be confused with the Legendre polynomials. The penetration factor and shift factor are
defined as [11]
P (E ) =

ka
2

F (ka) +G 2 (ka)

S (E ) = P (E ) F (ka)F (ka) +G (ka)G (ka)

(1.17)
(1.18)

where F and G are the regular and irregular Coulomb functions, F and G are their derivatives, k is the
wavenumber, and a is the channel radius, a somewhat arbitrary boundary chosen such that for r > a,
the scattering potential is approximately the Coulomb potential.
More practically, the R matrix theory can be used to calculate the collision matrix Ucc between channels c and c . In the Breit–Wigner approximation, this is
Ucc (E ) = e i (φc +φc ) δcc +

8

i Γc (E )Γc (E )
,
E R − E − i Γ(E )/2

(1.19)

where φc is the hard-sphere phase shift in this channel [11]. This can then be used to find the elastic
scattering cross section
dσ
1
=
d Ω (2I 1 + 1)(2I 2 + 1) M1 M2 M

(cM M )

1 M2

fC (Ω)δM1 M1 δM2 M2 + f c M 1M 2 (Ω)
1

2

2

(1.20)

with Coulomb scattering amplitude fC (Ω) and elastic scattering amplitude [11]
(c M M )

f c M 1M 2 (Ω) =i
1

2

π
k

Jπ I

I

2 + 1e i (σ

+σ )

× (I M 0|J M ) I 1 I 2 M 1 M 2 I M
× δc c δ I I δ

−UcJ πI

,c I

(I 1 I 2 M 1 M 2 |I M )
I

M M −M JM

(1.21)

Y M −M (Ω).

In these equations, I i for i = 1, 2 are the spins of the two colliding particles,

is the orbital angular mo-

mentum of the system, J is the total angular momentum of the many-body system, M is the projection
of the total angular momentum, and π is the total parity of the system. More detail about the R matrix
theory can be found in the review article by Descouvemont and Baye [11].

1.1.7 Isobaric analogue resonances
In some proton-adding or neutron-removing [22] reactions, isospin coupling between the target and the
projectile allows the excitation of higher-isospin states in the compound nucleus. For a target nucleus
beginning in a state with isospin T , these higher-isospin states have isospin T> ≡ T + 12 . Each T> state
is the isobaric analogue of states in other nuclei with the same mass number A but different isospin
projection T z ; therefore, a measurement of the properties of the resonance in the compound nucleus
allows the properties of the analogous states in the other isobars of the multiplet to be deduced.
These T> isobaric analogue resonances coexist with a potentially large number of T< ≡ T − 12 resonances corresponding to unbound states of the compound nucleus itself. If the density of the T< states
is large, the T> states may begin to mix weakly with T< states that have the same spin and parity [22].
This splits the T> resonance into several components, creating fine structure in the excitation function
as seen in Fig. 1.4. However, if the density of T< states is very high (or the resolution is lower), the individual components will not be resolvable, and the mixing will instead present itself as a broadening of
the T> resonance [22].

9

Figure 1.4: A high-resolution measurement of the excitation function from 40 Ar(p, p ). The bottom plot
shows a theoretical fit to the data near 2.45 MeV. Spin-parity assignments were made based on the shape
of the resonances, and it is apparent that the 1/2+ and 3/2− resonances are spit into several fine-structure
components. (Figure reproduced from [22, Fig. 19.13], which was adapted from [7].)
IAS
Er
δ
∆E c − δ

Sn

46

Ar + n

46

47

Ar + p

K

47

Ar

Figure 1.5: Calculating the resonance energy of the isobaric analogue state labeled “IAS”. In this figure,
S n is the neutron separation energy, δ is the difference between the neutron and proton masses, ∆E c is
the Coulomb energy of the last proton, and E r is the resonance energy. The ground state of 47 K is located
far below the bottom of the figure and is not shown. The energy separations are not drawn to scale.

10

1.1.8 Calculating the resonance energy
The expected resonance energies can be calculated by comparing the isobaric analogue state to similar
systems. Consider, for example, the state in 47 Ar shown on the right-hand side of Fig. 1.5. If we assume
that the strong nuclear force binds protons and neutrons equally, then the energy difference between the
state in 47 Ar and its analogue state in 47 K must be ∆E c − δ, where ∆E c is the Coulomb energy of the last
proton in 47 K, and δ = m n − m p is the mass difference between the proton and the neutron.
A similar comparison can be made between the bound 47 Ar nucleus and the unbound 46 Ar+n system;
in this case, the energy difference is simply the neutron separation energy S n . This unbound system has
an energy that is higher than the unbound 46 Ar + p system where the resonance occurs by an offset of δ.
Putting these pieces together, the energy of the proton emitted when the compound nucleus decays
from the isobaric analogue state is the difference in energy between the isobaric analogue state in 47 K
and the energy of the unbound 46 Ar + p system. This is
E r = (∆E c − δ) − S n + δ = ∆E c − S n .

(1.22)

Values of ∆E c can be determined empirically or by simulation; values of S n are typically found from mass
measurements and can be looked up in standard references.

1.2 Experimental motivation
The experiment described in this thesis served several purposes. From the perspective of nuclear physics,
the results are interesting because they can contribute to the understanding of the breakdown of the
magic numbers near the N = 28 shell closure. From a practical perspective, the experiment also served
to commission the Active-Target Time Projection Chamber (AT-TPC) with radioactive beams and to test
the analysis methods developed for AT-TPC data. These motivations are each described further in this
section.

1.2.1 The disappearing N = 28 shell gap
While the shell model and magic numbers describe stable and near-stable isotopes relatively well, the
shell gaps defined by the magic numbers begin to change in nuclei that lie further from stability. An

11

4000
800
3000

Ca
600

2

4

B(E2)up (e fm )

2000

+

E (2 1) (keV)

20

400

Mg

12

S

16

Ar
S
16
18

1000
12

Mg

Ar

200

18

Si

14

Si

14

0

0
18 20 22 24 26 28 30 32
Neutron number

Ca

20

18 20 22 24 26 28 30 32
Neutron number

Figure 1.6: Trends in E (2+
1 ) and B (E 2) near the N = 20 and N = 28 magic numbers. (Reproduced from
[49].)
example of this is predicted in the N = 28 isotones, where a reduction of the previously mentioned spinorbit splitting of the 1 f states could move the 1 f 5/2 level close enough in energy to the 2p 3/2 level to

γ

γ levels [18]. This would be indicainduce a strong electric quadrupole transition strength between the
tive of the onset of nuclear deformation, particularly of quadrupole nature [18], and the emergence of
collective behavior and excitations in these nuclei [18, 49].
Several
experiments have indicated that the N = 28 shell gap disappears as the number of protons
γ
decreases from 48 Ca. One piece of evidence cited by Sorlin and Porquet [49] is the evolution of the energy
of the first 2+ excited states E (2+
1 ) in these isotones as a function of Z . Magic nuclei that follow the shell
model would be expected to have large values of E (2+
1 ) due to the large shell gap just above their highestlying single-particle states. This is observed, for instance, in the N = 20 isotones, where E (2+
1 ) has a strong
peak in the doubly magic 40 Ca that decays very slowly through 38 Ar, 36 S, and 34 Si before disappearing in
32

Mg (see the left-hand plot in Fig. 1.6). Sorlin and Porquet [49] cite this as evidence of an N = 20 shell

closure that remains strong with decreasing Z until Z = 14. In the N = 28 isotones, on the other hand,
the strong peak in 48 Ca is already greatly attenuated at 44 S and completely gone by 42 Si. This indicates
that the N = 28 shell closure vanishes much more rapidly with decreasing Z [49].
A second piece of evidence for the departure from standard shell-model behavior is the reduced elec-

γ

+
tric quadrupole transition strength B (E 2; 0+
1 → 21 ) in these nuclei. This transition strength is expected

12

to be weak in spherical nuclei due to the large shell gap. As was the case with values of E (2+
1 ), this is
largely true for the N = 20 isotones, which all have a strong minimum for B (E 2) except for 32 Mg [49]. The
N = 28 isotones, on the other hand, have widely varying values of B (E 2), and the value for

46

Ar is not

even well-determined since different experimental techniques have yielded incompatible results [49].
This is shown in the right-hand plot of Fig. 1.6.
These measurements, along with β-decay and Coulomb excitation studies, have indicated that nuclear deformation develops smoothly as one removes protons from the spherical, doubly magic 48 Ca [49].
In particular, Sorlin and Porquet [49] note that prolate-spherical shape coexistence has been observed
in 44 S, 42 Si is known to be oblate, and 40 Mg is predicted to be prolate. This leaves 46 Ar in an interesting
shape transition region between the deformed 44 S and the spherical 48 Ca.
In addition to yielding information about nuclear shape in this region, studies of 46 Ar and 47 Ar can
also provide important information about single-particle structure near the N = 28 shell gap. Because
46

Ar has a closed neutron shell, the shell model predicts that the valence neutron in 47 Ar should occupy

the next available orbital, which is 2p 3/2 . Therefore, reactions populating the ground state and first excited state of 47 Ar should exhibit strong p-wave components, as revealed in their spectroscopic factors.
Additionally, measurements of the energy gap between the ground state of 47 Ar and its 1/2− first excited
state provide information about the magnitude of the spin-orbit splitting of the 2p harmonic oscillator
shell which, as mentioned previously, is expected to shrink as the N = 28 shell gap vanishes.
1.2.2 Previous studies of 47 Ar
The first [18] experiment to investigate the structure of
al. [16, 18]. They populated excited states in

47

47

Ar was performed in 2006 by Gaudefroy et

Ar via the d

46

Ar, 47 Ar p one-neutron transfer reaction

with a 10 MeV/u beam of 46 Ar generated by the SPIRAL facility at GANIL. The spectrum of 47 Ar (shown
in Fig. 1.7) was reconstructed from the energies and scattering angles of the outgoing protons, which
were measured using position-sensitive silicon telescopes at backward angles between 110° and 170°.
By comparison to shell model and distorted-wave Born approximation (DWBA) calculations, they were
able to make angular momentum and parity assignments to the ground state and first few excited states
(see Table 1.1). In addition, they assigned a neutron single-particle state of p 3/2 for the ground state and
tentatively assigned arrangements of p 1/2 and f 7/2 to the first two excited states at 1130 keV and 1740 keV,

13

(2J + 1)C 2 S

Jπ

E [keV]
0
1130(75)
1740(95)
2655(80)
3335(80)

3/2−
1/2−
7/2−
5/2−
5/2−

1
1
3
3,(4)
3,(4)

2.44(20)
1.62(12)
1.36(16)
1.32(18)
2.58(18)

Counts / 150 keV

Table 1.1: First few levels in 47 Ar identified by Gaudefroy et al. [16, Table 1].
140

b)

a)

100

120
100

50

80
0
c)

60

100

40
C+d

50

20
C
0
-4

-2

0

2

4

6

8

-4

-2

0

2

4

6

8

0

Excitation Energy (MeV)

Figure 1.7: Energy spectrum of d 46 Ar, 47 Ar p observed by Gaudefroy et al. [16]. Subplot (a) shows the
spectrum with background curves from carbon-induced reactions only (labeled “C”) and from carboninduced reactions and deuteron breakup (labeled “C + d”). Subplot (b) shows the background-subtracted
spectrum fit with 9 Gaussian curves. Subplot (c) shows the spectrum of events where protons were in
coincidence with the 47 Ar nucleus. (Reproduced from [16, Fig. 1].)
respectively [18]. Finally, they calculated spectroscopic factors for each measured state using a variety of
different optical potentials.
Using the measured energy levels, Gaudefroy et al. [16] calculated single-particle energies for the p f
shell states. They concluded that the spin-orbit splitting between the 2p levels in 47 Ar was reduced by
890(120) keV as compared to 49 Ca, and the splitting between the 1 f levels was reduced by 875(130) keV
(see Fig. 1.8). This produces a reduction in the N = 28 shell gap of 330(90) keV. This result is controversial,
however, and was disputed in a comment published by Signoracci and Brown [48]. They pointed out that
when computing the single-particle energies, Gaudefroy et al. [16] neglected to include hole states in
the weighted sum. Correcting this leads to a much smaller change in the p-state spin-orbit splitting of
10(130) keV [48]. In response to this comment, Gaudefroy et al. recalculated their values and published a
smaller value of 207 keV for the reduction in the spin-orbit splitting of the p states [17].

14

SPE (MeV)

-0 1

0

f5/2

-1.15

-2

-2.42

p1/2 -3.13

-3.55

-4

p3/2 -5.15

28

-6

28
-8.02

-8

-10

f7/2

-9.95
49

47

Ca

Ar

Figure 1.8: A comparison of neutron single-particle energies in 49 Ca and 47 Ar from Gaudefroy et al. [16].
The splitting of the 1p orbitals is reduced in 47 Ar compared to 49 Ca. (Reproduced from [16, Fig. 3].)
A second experiment performed by Bhattacharyya et al. [6] in 2008 measured the γ-ray spectrum
of

47,48

Ar produced by deep inelastic transfer reactions between a

238

U beam and a

48

Ca target. Their

results largely agreed with the previous measurement by Gaudefroy et al. [18], although as the reaction
mechanism was different, they populated a slightly different set of states. The levels identified in their
paper are shown in Fig. 1.9. Of particular interest is the new 5/2− state at 1234(4) keV which is very close
to the 1/2− first excited state. This was not observed by Gaudefroy et al. [16], and based on the resolution
of the peaks in Fig. 1.7, it would not have been possible for them to separate it from the 1/2− state.
Most recently, a thorough study of the levels in 47 Ar was presented by Gade et al. [15] in 2016. They
measured the γ-ray spectrum of 47 Ar at the NSCL using two different reactions. The first part of the experiment used the GRETINA array and the S800 spectrograph to measure the single-neutron pickup reaction
12

C

46

Ar, 47 Ar + γ X . The second part measured the reaction 9 Be

48

K, 47 Ar + γ X , or the removal of one

proton from the ground state of 48 K, using the SeGA array. These two measurements complemented each
other by populating different sets of excited states in the 47 Ar reaction product, leading to the two spectra shown in Fig. 1.10. Between the two measurements, the authors observed all bound states of the 47 Ar
nucleus that were known at the time of writing except for the 1/2− state at 1130 keV found by Gaudefroy

15

2.2
2.1

3/2 –
7/2 –

1.6
(5 / 2 – ) 1.3
(1 / 2 – ) 1.3

7/2 –
1/2 –
5/2 –

2190
516
1234

0

(7 / 2 – )

1200

2190

1234
1200

1747

1747

EXP

3/2 –

0

SM

3/2 –

47
18Ar29

Figure 1.9: Energy levels identified by Bhattacharyya et al. [6]. (Adapted from [6, Fig. 3].)
(a) 3438

5/2-

1693

2761
5/2-

1745
7/2515

47

1693

1017
3438

(b) 3438
5/2-

Ar from 46Ar+1n

47

1530

2188
3/2-

Ar from 48K-1p
2761
5/2-

1017

956

3438

1231
5/2- 1130
1/2- 2188
1745
1231

1745
7/2515

1530

2188
3/2-

956

1231
5/2- 1130
1/2- 2188
1745
1231
3/2-

3/2-

Figure 1.10: γ-ray spectrum measured by Gade et al. [15]. The results from each reaction are shown separately, and levels indicated in gray were not observed using a given reaction. (Reproduced from [15,
Fig. 3].)
et al. [18], which could not be populated due to orbital angular momentum constraints.

1.2.3 Resonant proton scattering in inverse kinematics
The resonant proton scattering method described in Section 1.1.7 was used very successfully in the past
in direct kinematics experiments. One particularly relevant example is the measurement of

40

Ar(p, p )

done by Scott et al. [47] in 1968. They bombarded a sealed nat. Ar gas target with a beam of protons from a
5.5 MeV Van de Graaff accelerator. The resulting excitation functions were measured to a high resolution
of approximately 3 keV, revealing the fine structure resonance components described in Section 1.1.7;

16

the gross structure of the resonances was recovered by averaging the data over approximately 25 keV.
Part of the reason these direct kinematics experiments were so successful was due to the characteristics of the proton beams they used. Being nothing more than 1 H, proton primary beams can be produced quite easily and with very high intensity as compared to radioactive secondary beams produced
by fragmentation. Furthermore, by virtue of being a primary beam, the incoming proton energies can be
controlled very precisely, producing a very small uncertainty in the reaction energy. The combination of
fine energy resolution and excellent statistics generally leads to impressive results, as seen in Scott et al.
[47].
A natural question, then, is whether this technique can be successfully extended to the modern
domain of inverse kinematics experiments which tend to feature a low-intensity beam with a slightly
broader energy profile. Some resonant proton scattering experiments have been performed in inverse
kinematics, but the technique has not yet been used with heavier radioactive beams in a device like the
AT-TPC. This experiment thus also serves as a test of the feasibility of this sort of measurement in this
detector.

1.2.4 Commissioning of the AT-TPC
In addition to the nuclear physics goals, the experiment described in this thesis was designed to commission the AT-TPC for use with radioactive beams. The AT-TPC is a newly built time projection chamber
with an active-target design that uses the detector gas as a scattering target. This allows reactions to occur
anywhere within the active volume of the detector, thereby increasing efficiency and allowing multiple
reaction energies to be measured at once. The AT-TPC is described in detail in Chapter 2.
Prior to this experiment, the AT-TPC was used to measure α-α elastic scattering at energies less than
4 MeV/u. The scattering angles measured in this experiment were reconstructed to a resolution of approximately 1°. Immediately prior to the 46 Ar experiment, 40 Ar(p, p ) elastic scattering was measured with
a stable 40 Ar beam in hopes of reproducing the results of Scott et al. [47], but an unforeseen issue with
gas contamination made this dataset more difficult to analyze. Therefore, this dataset will be processed
at a later time.

17

1.3 Summary
The remainder of this document describes an experiment that was performed to measure the reaction
46

Ar(p, p ) in inverse kinematics with the AT-TPC. As described above, this was done with the dual pur-

poses of commissioning the AT-TPC with radioactive beams and furthering the understanding of the
disappearance of the N = 28 shell closure. First, the AT-TPC and its instrumentation are described in detail in Chapter 2. The methods used to analyze AT-TPC data are then defined in Chapter 3. These analysis
methods are applied in both Chapter 4, which describes simulations of the experiment and the detector’s
efficiency, and Chapter 5, which presents the setup of the experiment, its data processing scheme, and
its results. Finally, a summary is given in Chapter 6 along with some notes about improvements currently
being made to the AT-TPC and its instrumentation.

18

CHAPTER 2
THE ACTIVE-TARGET TIME PROJECTION CHAMBER

As mentioned previously in Section 1.2.4, one of the goals of the experiment described in this document
was the commissioning of the Active-Target Time Projection Chamber (AT-TPC), a new detector that was
designed to measure reactions involving low-energy radioactive beams. This chapter begins with a brief
background on time projection chambers, while the remaining sections are devoted to a description of
the AT-TPC and the equipment used to instrument it during the experiment.

2.1 Background
2.1.1 Time projection chambers
The time projection chamber (TPC) is a type of gas-filled charged particle detector originally invented
by Nygren [42] in the 1970s. A TPC consists of a gas-filled vessel equipped with an anode and a cathode
capable of producing a moderate electric field. When a charged particle passes through the gas-filled
volume, it ionizes some of the atoms that compose the gas, leaving behind a trail of free electrons. The
electric field pushes these electrons toward the anode of the chamber, which is typically equipped with
some sort of position-sensitive amplification and readout device to count the number of electrons that
strike it. This process is illustrated in Fig. 2.1.
The position-sensitive readout plane gives two-dimensional information about the track; the third

Particle track

Anode

Cathode

E
Figure 2.1: Principle of operation of a TPC. A uniform electric field is created between the anode and the
cathode. This field transports the ionization electrons produced by a charged particle towards the anode,
where they are amplified and collected.

19

dimension is reconstructed from the timing of the recorded signals. This is possible because the electrons drift at a constant velocity in the detector [32], so the position along the drift direction is directly
proportional to the drift time. This is the origin of the “time projection” portion of the detector’s name.

2.1.2 Active targets
The TPC was originally designed as an evolution of bubble chambers for tracking particles in high-energy
experiments [32, 42]. To observe nuclear reactions inside a TPC, a sample of target material must be installed within the detector. This sample is traditionally a thin foil of target material or a small amount of
target material evaporated onto a thin substrate, but the target can also be the TPC gas itself, creating
a so-called active target [4]. This design has been used successfully in several recently constructed detectors for nuclear physics applications, some examples of which are described by Beceiro-Novo et al.
[4].
The active target design provides a number of advantages over a solid target. The main benefit is that
it allows the detection of reactions at low energy with a solid angle coverage close to 4π and without
compromising the target thickness, as opposed to a passive target where the thickness has to be reduced
enough for the reaction products to escape. Because reactions can occur anywhere within the detector’s
active volume, they can be measured over a continuous, broad range of energies as the beam slows down
within the gas volume. Since the reaction products emerge from within the tracking medium itself, very
low energies can be measured, and the vertex of the reaction can be reconstructed on an event-by-event
basis, allowing excitation functions to be measured with a single beam energy.
The active target design imposes some constraints on the choice of gas for the TPC. Not only does the
gas have to provide good electron amplification and transport, but it also must contain the target nucleus
of interest and, ideally, few other nuclei. For example, the proton scattering experiment described in this
thesis required a gas containing hydrogen nuclei. The best gas choice would therefore be pure hydrogen
as it has no other nuclei for the beam to scatter off of. However, the electron amplification properties of
hydrogen make it hard to work with. Instead, we chose to use isobutane (C4 H10 ). While isobutane provides better electron amplification, it also requires us to filter out scattering off of carbon nuclei during
the analysis.

20

AsAd boards

Sensor Plane
& Micromegas
Cathode

High-voltage
feedthrough

Beam
entrance

Field cage
Beam duct

Figure 2.2: A schematic view of the AT-TPC. The outer shielding volume was made transparent in this image to make the details of the inner volume more visible. Beam enters the detector through the beam duct
at the right-hand side of the image and moves toward the sensor plane on the left. Some components of
the GET electronics are shown mounted on the downstream end of the detector (see Section 2.4).

2.2 Structure of the AT-TPC
The AT-TPC is a newly constructed time projection chamber that uses the active target design described
above. Figure 2.2 shows a cutaway view of the detector and all its components. The structure of the TPC
is described briefly in this section. A more thorough description of the detector’s design can be found in
an article published by Suzuki et al. [50] describing a half-scale prototype of the AT-TPC.

2.2.1 Inner and outer gas volumes
The active volume of the AT-TPC is formed by an open-ended cylinder of epoxy-coated fiberglass with
length 1 m and radius 28 cm. The upstream end of the volume is sealed with a stainless steel cathode,
which has at its center a thin foil window through which the beam enters the detector. The downstream
end is closed by an aluminum flange which supports the sensor plane that serves as the anode. This
inner volume is filled with the gas that serves as a target.
Surrounding the inner volume is a larger, concentric shielding volume contained by a cylindrical
aluminum pressure vessel. The purpose of the shielding volume is to isolate the high electric potential of
the cathode from the environment and to prevent arcing. Therefore, the shielding volume is filled with

21

Figure 2.3: A photograph of the inner volume wall with the rings of the inner field cage installed. The
outer field cage rings have a larger diameter, but are otherwise identical and are installed on the outer
surface of the cylinder.
an inert gas with a high dielectric constant, like nitrogen. The pressure of the shielding volume is kept
slightly lower than the pressure in the active volume to ensure that if there were a small leak between
the two volumes, the inner volume would not be contaminated with a foreign gas that could affect the
electron drift velocity.

2.2.2 Electric and magnetic fields
The active volume contains a uniform electric field which is produced by applying a potential difference
between the cathode at its upstream end and the anode at its downstream end. To ensure that the electric
field is uniform, the walls of the inner volume are surrounded inside and out by a field cage consisting of
concentric ring-shaped electrodes connected by a resistor chain (see Fig. 2.3). The resistor chain gradually steps down the voltage between each ring, ensuring electric field uniformity throughout the volume.
In addition to the electric field, particles in the AT-TPC experience a magnetic field from the largebore solenoidal magnet that the detector is mounted inside of. This magnet was originally designed for

22

use in a medical magnetic resonance imaging (MRI) machine, and it is capable of producing a field of up
to around 2 T. The AT-TPC is mounted on rails in the center of the magnet. The longitudinal magnetic
field bends the trajectories of the emitted charged particles in order to better determine their energies.
Another benefit is the ability to track particles over longer paths, and for those that stop in the gas volume,
measure their total range.

2.2.3 Sensor plane
Ionization electrons produced inside the AT-TPC are read out by a sensor plane mounted on the anode
end of the detector. The sensor plane (or “pad plane”) consists of a circular printed circuit board of radius
27.5 cm covered with 10 240 triangular gold-plated conductive pads. The pads are arranged in a hexagonal inner region of small pads with height 0.5 cm surrounded by an outer region of large pads with height
1.0 cm. This layout is shown in Fig. 2.4.
The triangular pad shape was chosen to maximize the spatial resolution of the detector. When a
charged particle track crosses a series of adjacent triangular pads, the amount of charge deposited on
each pad is staggered. This staggering pattern changes significantly with even a small change in the orientation of a track, leading to a greater sensitivity and an improved angular resolution.
Electron amplification is provided by a Micromegas [19] device installed on the sensor plane. The Micromegas (from “micro-mesh gaseous structure” [27]) consists of a very fine conductive mesh supported
approximately 100 µm above the sensor plane by insulating posts. The mesh is biased with respect to the
electrodes to create a relatively large electric field between the mesh and the pads. This effectively divides
the detector into two regions: a large drift region above the mesh and a small multiplication region below
the mesh. When a charged particle passes through the detector, it ionizes the gas in the drift region. The
relatively low electric field in the drift region then transports the electrons toward the mesh. Once the
electrons pass through the mesh and enter the multiplication region, the much larger electric field there
causes electron avalanches to form, effectively amplifying the single-electron signals to something that
can be measured by the electronics. This process is illustrated in Fig. 2.5.

23

Figure 2.4: Layout of the pad plane. The inset shows a closer view of one corner of the hexagonal inner region. This region of half-scale pads provides finer resolution near the reaction vertex, which will
generally occur near the central axis of the detector.

Cathode

Particle
E

5

10 V/m
e

−

e−

e−
Mesh

E ∼ 106 V/m

Pads
Anode

Figure 2.5: Principle of operation of a Micromegas. The electric field magnitudes shown are nominal: the
field in the region above the mesh may vary from roughly 103 V/m to 105 V/m.

24

Figure 2.6: Illustration of pileup separation due to the tilt angle. The leftmost track is elastic α-α scattering, and the rightmost track is another α particle that entered the detector during the event. Without the
tilt angle, the two beam tracks would be collinear, making them more difficult to separate.

2.3 Tilting the detector
The AT-TPC is normally coaxially centered within the bore of the magnet, but the detector’s support
carriage has a small jack that allows the anode end of the detector to be raised up to approximately 6°
with respect to the axis. Tilting the detector can be advantageous for two reasons. First, the tilt angle
projects tracks at very forward angles onto more than one pad, improving measurements at these small
scattering angles. Second, it allows some separation of pileup events since earlier beam tracks drift away
from the cathode before later ones arrive (see Fig. 2.6).
While tilting the detector solves some problems, it also creates a few complications for the analysis.
When tilted, the electric and magnetic fields in the TPC are no longer parallel, but are instead separated
by a small angle. This causes the drift electrons to be deflected by a small amount in a direction transverse
to the electric field. This has to be corrected for during the track reconstruction process, which will be
described in Section 3.2.

25

MuTAnT

µTCA Backplane

CoBos

AsAd

10 Gb/s Fiber-optic link

CoBos

Gigabit Ethernet Switch

Control &
Configuration

ADC

DAQ Cluster

AGET AGET AGET AGET

Long-term storage
Isolation circuit
Sensor plane

Figure 2.7: A schematic view of the GET electronics system. Only one AsAd is shown for clarity, but 40
AsAds are used to instrument the AT-TPC.

2.4 Electronics
The AT-TPC is instrumented with digital electronics developed by the Generic Electronics for TPCs (GET)
collaboration [44]. This equipment provides a fully digital data acquisition system capable of digitizing
and recording the full trace for each of the 10 240 channels in the detector.
The electronics hardware is divided into a hierarchy of several modules (see Fig. 2.7). At the lowest
level of the hierarchy is a custom application-specific integrated circuit (ASIC) called the AGET (ASIC for
GET). This chip samples and shapes the signals and compares them to a threshold to generate a channellevel trigger. The AGET amplifies the incoming signal with a variable-gain charge-sensitive preamplifier
and performs pole-zero correction. It then stores samples of the analog signal in a switched capacitor
array (SCA) which acts as a circular buffer [1, 3]. Each AGET can read out 64 physics channels from the
detector in addition to 4 fixed-pattern noise channels. The fixed-pattern noise channels are structurally
identical to the physics channels, but without inputs [1, 3]. This allows the data to be corrected for lowfrequency electronic noise.
The AGETs are mounted in groups of four on AsAd (ASIC Support and Analog to Digital conversion)
boards. In addition to the four AGETs, each AsAd board houses a four-channel, 12-bit ADC. When a trigger is issued, the ADC digitizes the sample outputs from each AGET chip and transmits them via a serial

26

link [3]. Between triggers, the ADC digitizes and transmits the multiplicity signal from each AGET. The
input end of each AsAd board is attached to an isolation circuit, and this assembly is then mounted on
the downstream end of the TPC and connected directly to the sensor plane (see Fig. 2.8). This design
was chosen to minimize the distance that analog signals must travel before being digitized, reducing the
capacitance and potential for noise in the data.
The top level of the GET electronics hierarchy is the CoBo (Concentration Board). Each of the ATTPC’s 10 CoBos is connected to four AsAd boards. When a trigger is issued, the CoBo collects the data
from these boards, applies an event time stamp, and builds the event [3]. It then sends the event over a
10 Gb/s fiber-optic link to a network switch to be distributed to a cluster of computers for storage.
To keep the CoBos synchronized and generate a global trigger, an additional board called the MuTAnT (Multiplicity, Trigger, And Time) distributes a global time stamp and manages clock synchronization across the system [3]. The MuTAnT board is designed to collect all of the running multiplicities
and hit patterns from the CoBos and combine them in various ways to generate a global trigger. The
most straightforward method consists of simply summing the multiplicities to generate a trigger when a
global threshold is reached. This is referred to as a Level 1 trigger. Other more complex trigger decisions
can be used that involve matching a predefined pattern of hit channels (a Level 2 trigger). The choice
of which channels contribute to the global multiplicity trigger is set independently for each CoBo; this
process will be described in Section 5.2.3.

2.5 Data acquisition
Due to the very large number of channels and high data rate, the AT-TPC requires custom software for
data acquisition. The AT-TPC DAQ is a multi-level, distributed system with a controlling master node
overseeing a set of 10 worker nodes, one per CoBo. A schematic of the system is shown in Fig. 2.9.

2.5.1 Structure of the DAQ system
Since the DAQ computers must be located next to the AT-TPC in the experimental vault, the DAQ system
was designed as a web application that can be opened in a browser on a client computer elsewhere in
the laboratory. The core of the software is written in Python 3 using the popular Django web application

27

Figure 2.8: A photograph of the fully instrumented AT-TPC mounted inside its solenoid magnet. The
40 AsAds are mounted on the downstream end flange surrounded by copper shielding. Individual cables
connect each AsAd to its controlling CoBo, and a separate set of cables connects each to its power supply.
framework and the Celery asynchronous task queue. The site is supported by a PostgreSQL database that
stores the internal configuration of the system.
The processes for the web servers controlling the DAQ interface are run inside a set of Docker containers. Docker is a containerization platform that allows programs and their associated dependencies
to execute independently and without knowing about the rest of the system [12]. This is possible through
the use of containers, which are functionally similar to virtual machines but with the key distinction that
a container does not contain a full guest operating system. Instead, containers use the Linux kernel of
the underlying host’s operating system while still providing filesystem isolation and an independent set
of system libraries to the contained code. This saves disk space compared to a full virtual machine, but
more importantly, code inside a container runs at native speed as it does not have to work through a
hypervisor layer. Additionally, since each container has its own full set of system libraries, containerized
applications are platform-independent and can run on any operating system that uses the Linux kernel
(support is also available for Windows and macOS through a lightweight Linux virtual machine).
In addition to the web application, there are ten separate worker nodes for the DAQ, one for each
CoBo. Each of these runs on its own computer. The main reason for running so many independent systems is throughput: with a theoretical maximum1 data rate of approximately 21 MB/event, the AT-TPC
1 Experimental data rates are much lower since we do not read out channels that were not above threshold.

28

Master node

CoBo config files

Docker
Controller

Message
Broker

Async
tasks

Client
Web browser

Web server
Database

Worker node

Worker node

Electronics control

Electronics control

Data router

Data router

Storage

Storage

CoBo

CoBo

…

Figure 2.9: Structure of the DAQ system. Only two worker nodes are shown, but there are ten in the actual
system.

29

could easily saturate the throughput of a single mechanical hard drive at an event rate as low as 10 Hz.
Using a parallel set of 10 computers allows the rate to be 10 times higher.
Each worker node runs a set of two programs provided by the Institut de recherche sur les lois fondamentales de l’Univers (IRFU) in Saclay, France. The first program, the Electronics Control Core (ECC)
server, is responsible for controlling the CoBo. It can configure the CoBo and its connected AsAd boards
using a configuration file and start and stop data taking. The second program, the data router, receives
the data stream from the CoBo and records it to a file on the local hard disk. At the end of the experiment,
the data files from the ten computers are collected and merged together.

30

CHAPTER 3
ANALYZING AT-TPC DATA

Extracting physics results from the data recorded by the AT-TPC requires a complex, multi-step analysis
process. This process involves correcting the electronics signals for baseline fluctuations, reconstructing
three-dimensional tracks, separating the tracks from noise, and fitting the tracks using a Monte Carlo
algorithm. This chapter contains a general description of the steps of the analysis process and the reasoning behind its design. Specifics of the analysis performed on the 46 Ar data are given later in Chapter 5.

3.1 Baseline correction
The signals recorded by the AT-TPC in the 46 Ar(p, p ) experiment feature large-scale baseline fluctuations
corresponding to the time that the beam particle was present in the chamber. This happened because
the 46 Ar nucleus, which is much more ionizing than the protons that the AT-TPC was tuned to measure,
deposited enough charge onto the pads it activated on the sensor plane to induce a negative signal in all
other pads via capacitive coupling through the mesh. Fortunately, this effect is easy to correct in software.
The correction is done using a Fourier transform and a low-pass filter.
One definition of the Fourier transform of a function f (t ) is1
fˆ(ν) = F [ f (t )] =

∞
−∞

f (t )e −2πi νt d t .

(3.1)

To correct for the baseline fluctuation, we first calculate the Fourier transform fˆ(ν) of the signal. Then,
we multiply fˆ(ν) by sinc(ν/a) where
sinc x ≡

sin(πx)
,
πx

(3.2)

and a is a constant scaling factor. Finally, we find the inverse Fourier transform of this product. This gives
an approximation for the baseline of the signal which can then be subtracted away (see Fig. 3.1).
The sinc function was chosen in order to take advantage of the Convolution Theorem. This theorem
states that
F −1 [ fˆ(ν)gˆ (ν)] = ( f ∗ g )(t ) =

∞
−∞

f (τ)g (t − τ)d τ,

(3.3)

1 This is perhaps not the most common definition in physics, but it is the definition used by the Fast Fourier Transform
library used in the analysis [25].

31

Amplitude [ADC bins]
Amplitude [ADC bins]

500
475
450
425
400
375

Raw signal
Baseline

100

Corrected signal

75
50
25
0
0

64

128

192

256

320

ADC time bucket

384

448

512

Figure 3.1: Subtracting the baseline from a single signal. The top plot shows the raw signal and the calculated baseline. The bottom plot shows the difference between the two.
where ( f ∗ g )(t ) denotes the convolution of f (t ) and g (t ) [31]. Therefore multiplying the transformed
signal by sinc(ν/a) and then inverting the transformation is equivalent to convolving the original signal
with the inverse transformation of sinc(ν/a). Conveniently, that inverse is

1
1


1, − < t <
2
2
F −1 [sinc(ν)] = rect(t ) ≡


0, elsewhere,

(3.4)

a rectangular window function. This whole procedure is therefore equivalent to convolution with a rectangular window, or a moving average. The Fourier transform was used instead of a simple moving average due to the high speed of modern Fast Fourier Transform (FFT) libraries.

3.2 Track reconstruction
After baseline correction, the next step of the analysis process is to reconstruct three-dimensional tracks
from the signals. This entails converting the time-domain signals to a set of discrete points (x i , y i , z i ) in
three-dimensional space with associated peak amplitudes a i . The simplest way to do this is to assume

32

y
τ

v

B
E

τ

B
E

x

w

z

u

τ

τ

Beam

Beam

(a) Detector coordinate system

(b) Beam coordinate system

Figure 3.2: Coordinate systems used in the analysis. After calibration, the recorded tracks are in the detector coordinate system (a). This system has its origin in the center of the sensor plane at z = 0 with the
z axis pointing upstream toward the beam entrance window. The beam coordinate system (b) is related
to the detector coordinate system by a rotation through angle τ about the beam entrance window, where
τ is the detector tilt angle. In this system, the beam moves along the −w axis. Note that both of these
coordinate systems are left-handed.
that each pad is activated only once in a particular event. Then the raw (x, y) coordinates for a particular
channel are assumed to be the centroid of the corresponding triangular pad, and the uncalibrated z
position is the location of the center of gravity of the peak in the ADC signal. The amplitude of the signal
is calculated by subtracting an averaged baseline from the maximum of the peak. This averaged baseline
is the mean of the signal in time buckets N − 20 through N − 15, where N is the location of the maximum
of the peak.
At the end of this process, we have a set of 4-tuples (x i , y i , z i , a i ) that describe the position and energy
loss of a tracked particle in the detector coordinate system (Fig. 3.2a) at a given moment in time; however,
the units of these quantities are not all the same. Since the x and y positions are calculated from the pad
plane, they are in distance units, but the z position is still an ADC time bucket index, which is a time unit.
A calibration is therefore required to make all of these values compatible with each other.
The calibration is simplest when the electric and magnetic fields are parallel or when there is no
magnetic field. In this case, the primary electrons produced in the TPC drift along the electric field lines,
which point along the z axis in the detector. Therefore the track’s projection onto the pad plane is orthogonal, which implies that the x and y coordinates recovered from this projection equal the true coordinates of the particle track, and the only thing that needs to be calibrated is the z position. Since the
uncalibrated z value is in time units, it can be converted to distance units using the electron drift velocity

33

v d and the ADC clock frequency ν:
z=

v d z0
.
ν

(3.5)

The calibration is more complex when the magnetic field is nonzero and the detector is tilted. In this
case, the electric and magnetic fields are neither parallel nor perpendicular, so the Lorentz force experienced by the drifting electrons will have components in the transverse x and y directions in addition to
the primary z component. This causes the electrons’ trajectories to be deflected in the x and y directions,
which means that all three coordinates need to be calibrated. This can be done by using a vectorial drift
velocity instead of the scalar one used above.
Following the derivation presented by Lohse and Witzeling [32], the motion of an electron with mass
m and charge e can be modeled with a Langevin equation of the form
m

m
dv
= e(E + v × B) − v,
dt
τ

(3.6)

where E and B are the electric and magnetic field vectors, v is the electron’s velocity, and τ is the mean
time between collisions. This equation admits a steady-state solution
vD =

µE
ˆ .
ˆ ×B
ˆ + ω2 τ2 E
ˆ ·B
ˆ B
ˆ + ωτ E
E
1 + ω2 τ2

(3.7)

Here, vD is the drift velocity vector, ω = eB /m is the cyclotron frequency, τ = mv D /eE is the mean time
between collisions, and µ = eτ/m is the electron mobility in the gas. In the detector coordinate system
ˆ y,
ˆ z)
ˆ shown in Fig. 3.2a, the electric and magnetic fields are
(x,
E = E zˆ
B = B sin(θt ) yˆ + cos(θt )zˆ
for a tilt angle θt , so Eq. (3.7) can be simplified to yield the following components:
µE
(ωτ sin θt )
1 + ω2 τ2
µE
vy =
ω2 τ2 sin θt cos θt
1 + ω2 τ2
µE
vz =
1 + ω2 τ2 cos2 θt .
1 + ω2 τ2
vx =

34

(3.8)
(3.9)
(3.10)

Since we assumed a constant drift velocity magnitude, the amount of deflection in the x and y directions
is proportional to the drift time. Therefore, the data can be calibrated with the equations
x = x0 −
y = y0 −
z=

v x z0
ν
v y z0
ν

v z z0
,
ν

(3.11)
(3.12)
(3.13)

where v x , v y , and v z are the three components of the drift velocity vector given in Eqs. (3.8) to (3.10).
The final step of track reconstruction is to transform the data from the detector coordinate system
x y z of Fig. 3.2a to the beam coordinate system uv w of Fig. 3.2b. The beam is parallel to the w axis in the
latter frame, which simplifies the fitting algorithms. This transformation is given by
  

  
0
 u  1
  
  
 v  = 0 cos θt
  
  
w
0 − sin θt

0

0

 x  

  

  
 y  + (1 m) tan θt  .
sin θt 

  

  
cos θt z
0

(3.14)

The extra factor added to the v component corrects the vertical offset introduced by rotating about the
x axis at z = 0 rather than at the beam origin window.

3.3 Noise removal
The reconstructed tracks generally contain some noise points and structures that need to be removed to
prevent fits from diverging. There are generally two types of noise present in the events: uncorrelated,
random points from channels that were triggered by electronic noise, and correlated structures of points
created by cross-talk in the electronics. An example event with both types of noise is shown in Fig. 3.3.
The noise is removed by two algorithms that are run in series on each event. First, correlated noise is removed using a method based on the Hough transform [13], and then uncorrelated noise is removed using
a nearest-neighbor algorithm. This section begins with a description of the Hough transform algorithm
for straight lines and for circles, and then the noise removal algorithms are defined.

3.3.1 Hough transform for lines
The Hough transform is a feature recognition algorithm commonly used in computer vision applications
to detect straight lines and simple curves in images [13]. It works by transforming image points into a
35

200
150

v [mm]

100
50
0
50
100
150
100

0

100

u [mm]

200

0

500

1000

w [mm]

1500

2000

Figure 3.3: Event 99-30 from the 46 Ar data set, a particularly noisy event, before cleaning. There are some
correlated noise structures present at the bottom of the uv projection and the left of the w v projection,
and there is uncorrelated noise throughout the event.
parameter space where the desired image features map onto single points. The features can then be
identified by using a peak-finding algorithm in the parameter space.
One way to parameterize a line is in terms of the normal line drawn from it to the origin. This parameterization is given by [13]
R = x cos θ + y sin θ,

(3.15)

where R is the length of the normal line, and θ is the angle between the normal line and the x axis. This
is illustrated in Fig. 3.4. This parameterization implies that each point in the (x, y) system (the coordinate
space) maps onto a sinusoidal curve R(θ) in the (θ, R) system (the Hough space). Similarly, each line in
the coordinate space has a unique parameterization (θ, R), mapping it onto a single point in the Hough
space. Therefore, sets of collinear points in the coordinate space can be found by looking for intersecting
curves in the Hough space, as shown in Fig. 3.5.
This process is described more concretely in Algorithm 1. The algorithm computes a discrete set of
points (θi , R i ) in the Hough space corresponding to the sinusoidal curve generated from each data point.
Each of these points is then mapped to the corresponding bin in a two-dimensional accumulator array,
which is incremented by 1. After processing each point, lines in the original data are found by looking for
peaks in the accumulator array. A large number of collinear points creates a set of sinusoids that intersect
at one point, causing the corresponding accumulator bin to hold a large value.

36

y

R
θ

x

Figure 3.4: Parameterizing a straight line in two dimensions using a normal line drawn to the origin.

y

R

5
4

BC

4
BC

AB

2
B

3

A

A
B

2

0.5

AB

1.0

1.5

2.0

2.5

3.0

θ

-2

1

C

C
1

2

3

4

5

x

-4

(a) Coordinate space

(b) Hough space

Figure 3.5: Transforming points into the Hough space. The left-hand plot shows three points A, B , and C
in the normal coordinate space with lines drawn between them. These points map to sinusoidal curves
in the Hough space, as shown in the right-hand plot. Lines in the normal system become points in the
Hough space, so a line that passes through two points in the normal space is the intersection of two
curves in the Hough space.

37

Algorithm 1 Hough space binning method
Given: an Nbins × Nbins accumulator array A, a set of (x, y) data points D, and a maximum bin value M ,
function H OUGH B INS(A, Nbins , D, M )
for i ∈ {0, 1, . . . , Nbins − 1} do
i is the θ bin index
θ ← i π/Nbins
for (x, y) ∈ D do
R ← x cos θ + y sin θ
j ← (R + M )Nbins /(2M )
j is the R bin index
Ai j ← Ai j + 1
Increment the Hough space accumulator array
end for
end for
return A
end function
y

R
θ

x

Figure 3.6: Finding the center of a circle using the Hough transform.
3.3.2 Hough transform for circles
Although the process described in Section 3.3.1 is specific to finding lines, a similar process can be used
to find circles in an image. This procedure was originally described by Illingworth and Kittler [23], but
the formulation shown in this section was taken from Heinze [21].
Consider a set of points lying along a circle. If a line is drawn between any two of the points, a second
line can be drawn perpendicular to the first line and passing through its midpoint. This second line will
always pass through the center of the circle. Therefore, given a set of points, we can find the center of
curvature by finding the point where all of these perpendicular bisectors intersect (see Fig. 3.6). We can
do this using the Hough transform.

38

Given two points on the circle (x 0 , y 0 ) and (x 1 , y 1 ), the midpoint between the two points is
xm =

x0 + x1
;
2

ym =

y0 + y1
.
2

(3.16)

The slope of the perpendicular bisector to the line between these two points is
m=

x1 − x0
,
y0 − y1

(3.17)

and the y-intercept is
b = y m − mx m .

(3.18)

These can be plugged into the formula for a line y = mx + b and simplified to give
y=

(y 2 − y 02 ) + (x 12 − x 02 )
x1 − x0
x+ 1
.
y0 − y1
2(y 1 − y 0 )

(3.19)

This describes the set of possible locations of the center of the circle. To make the calculation easier, we
can transform this to polar coordinates:
x = R cos(θ)

(3.20)

y = R sin(θ).

(3.21)

Then, we can rewrite Eq. (3.19) as
R=

(x 12 − x 02 ) + (y 12 − y 02 )

2[(x 1 − x 0 ) cos θ + (y 1 − y 0 ) sin θ]

.

(3.22)

This equation can then replace the line R ← x cos θ + y sin θ in Algorithm 1 to make a Hough transform
function that finds the center of a circle of points.

3.3.3 Removing correlated noise
The two Hough transform algorithms described above are used to remove correlated noise points from
the recorded events. This process proceeds as follows:
1. Find the coordinates (u c , v c ) of the center of curvature of the event using the Hough transform for
circles (Section 3.3.2).

39

2. Calculate the arc length coordinate for each data point as a function of w. The arc length coordinate for each point is defined as
s = ρφ

(3.23)

where
ρ=

(u − u c )2 + (v − v c )2

and φ = tan−1

v − vc
.
u − uc

(3.24)

This parameterization is used because the arc length swept by the particle increases monotonically
as a function of w due to the constant polar angle θ of the particle’s momentum vector.
3. Perform the Hough transform for lines (Section 3.3.1) on the ρφ vs. w projection of the data.
4. Extract the top five peaks from the Hough space and average their angle values together to find the
angle of maximum activation in the Hough space. (This takes advantage of the fact that the lines
we are searching for in the ρφ vs. w space should be parallel.) Take an angular slice of the Hough
space around this maximum value and look for peaks in the 1D projection of this slice onto the R
axis. These peaks correspond to individual lines in the original space.
5. Determine which of the lines each data point is closest to. If a point is more than some threshold
distance from any line, it is discarded as noise.
6. Finally, check how many points are assigned to each line. If any line has less than some threshold
number of points, the points assigned to that line are discarded as noise.
The results of this process for one event are shown in Fig. 3.7.

3.3.4 Removing uncorrelated noise
Random, uncorrelated noise points that are left after the Hough transform–based cleaning are removed
by counting the number of neighbors each point has and eliminating points with too few neighbors
inside a given radius. We use a naïve algorithm to find the neighbor counts which simply calculates the
pairwise distance between each unique set of points in the data set and increments a counter if the
distance is small enough. This is shown in detail in Algorithm 2. This method is probably not the most
efficient choice, but it runs fast enough for our purposes.

40

300

R [mm rad]

200
100
0
100
200
300
0

250

500

750

1000 1250 1500 1750 2000

w [mm]

(a) Classification of data points

500

Hough space R bins

400
300
200
100
0

100

200

300

Hough space bins

400

(b) Hough space with angular slice

Figure 3.7: Cleaning an event with the Hough transforms. Plot (a) shows the data from Fig. 3.3 in the ρφ vs
w space with the found lines. The colored points are classified as data, while the empty points are noise.
Plot (b) shows the Hough space calculated for this data. The maximum angular slice is shown between
the two white lines, and the 1D projection of this slice is shown on the right with the peak locations
marked.

41

200
150

v [mm]

100
50
0
50
100
150
100

0

100

u [mm]

200

0

500

1000

w [mm]

1500

2000

Figure 3.8: Final results of the noise removal process on the event shown in Fig. 3.3. The blue points
were kept, and the semi-transparent red points were removed as noise. Some noise is still present after
cleaning, but much of it was removed.
After running both noise removal procedures, much of the noise in the event is removed (see Fig. 3.8).
The remaining noise is often focused in the first or last few time buckets of the event, and can sometimes
be removed with a simple cut. Regardless, the level of noise remaining after this process seems to be
small enough that it does not affect the Monte Carlo fits too much.
Algorithm 2 Counting nearest neighbors inside a radius
Given data X ∈ RN ×M , with number of data points N and number of dimensions M , and maximum
neighborhood radius r max ,
function N EIGHBOR C OUNT(X , r max )
let A ∈ ZN contain the neighbor counts initialized to −1 to negate self-counting
for i ∈ {0, 1, . . . , N − 1} do
for j ∈ {0, 1, . . . , i − 1} do
r2 ←

M −1

(X i k − X j k )2

k=0
2
if r 2 < r max

then
Ai ← Ai + 1
end if
end for
end for
return A
end function

42

3.4 Modeling tracks and fits
The result of the noise removal process is a set of data points that describe the track as a series of threedimensional locations with associated signal amplitudes. A model must be fit to this track in order to
extract physical observables like the scattering angle or the energy of the beam particle at the vertex.
This fit is performed using a Monte Carlo–based optimization algorithm, which will be discussed at the
end of this section.
Note that in the remainder of this chapter, the coordinates of the particle will be referred to using the
variables x, y, and z rather than u, v, and w. Regardless of the notation, this should be understood as a
position in the beam coordinate system shown in Fig. 3.2b.

3.4.1 Track representation and generation
Particle tracks in the AT-TPC can be uniquely described by a set of six parameters that define the initial state of the tracked particle. These parameters are the particle’s initial position (x 0 , y 0 , z 0 ), its initial
kinetic energy E 0 , and the azimuthal angle φ0 and polar angle θ0 defining the orientation of its initial
momentum vector in spherical coordinates. The azimuthal angle is measured with respect to the x axis
of the detector, and the polar angle is measured with respect to the z axis.
Given this set of parameters, a particle track can be modeled by simulating the motion of the particle
in the TPC under ideal conditions. The simulation is iterative, beginning with the set of initial coordinates described above. At each time step, it updates the position, momentum, and energy of the particle
after applying the Lorentz force. The updated energy is used to calculate the amount of energy lost to
interactions with gas atoms or molecules during that time step, and this energy loss is subtracted from
the updated energy. Finally, if the particle’s energy has reached zero, the simulation is stopped.
Ideally, it would be more convenient to model tracks with a mathematical function rather than a simulation; however, the nature of the AT-TPC makes this impossible. The particle tracks that the detector
produces are essentially helical, but the radius of curvature of the helices decreases in a complex way that
depends on the energy loss function, a Bragg curve. This prevents us from using an analytic function to
describe the tracks.

43

y

σ

x
Figure 3.9: Template used to simulate diffusion. The charge q from the initial point is distributed among
the diffused points. The central point gets 0.4q units of charge, and the outer points each get 0.6q/8 units
of charge. This is based on a two-dimensional normal distribution.
3.4.2 Simulated diffusion and track projection
The simulated track represents the path the particle followed through the inner volume of the TPC. To
perform the fit, we also need to calculate the distribution of charge deposition on the pad plane, or the
hit pattern. This is done in two phases.
First, we simulate electron diffusion. As electrons drift toward the end of the detector, the electron
cloud spreads out laterally and attains a final width of
2D t ,

σ=

(3.25)

where t is the total drift time (which is directly proportional to z), and D is a constant that is determined
empirically [27]. Diffusion is estimated by splitting each simulated track point into nine diffused points
based on the template shown in Fig. 3.9.
Each of the diffused points is then projected onto the readout plane at z = 0 using the inverses of
Eqs. (3.11) to (3.13), which were used to calibrate the track. A fine-grained two-dimensional lookup table
is used to map simulated points to channels in the electronics. The charge associated with each point is
then used to simulate a pulse in the appropriate channel, and this pulse is accumulated into an eventwide signal for that channel. The pulses are generated with the equation
f (q, t ) = A(q)

t
s

3

e −3(t /s) sin

t
s

(3.26)

with shaping time s and amplitude
A(q) =

Gµ q
G elec. NADCC norm
44

,

(3.27)

where G µ is the gain of the Micromegas, G elec. is the gain of the GET electronics, q is the amount of charge
deposited by this point, NADC is the total number of ADC bins, and C norm is a normalization factor given
by
C norm = max

f (q, t )
= max
A(q)

t
s

3

e −3(t /s) sin

t
s

= 0.044.

(3.28)

After all of the points are simulated, the peak amplitude is found in each channel to directly compare to
the experimental data.

3.4.3 The objective function
To compare data to the model, we need an objective function that describes how well the modeled track
fits the data. The objective function contains three contributions: one from the particle’s position, one
from its energy loss, and one that helps constrain its vertex location.
The position component of the objective function is a simple least-squares comparison between the
modeled track and the experimentally measured track. The modeled track’s x and y components are
˜
˜
linearly interpolated as a function of z to give continuous functions x(z)
and y(z).
These functions are
then evaluated at each of the experimental track’s z positions and compared to the experimental x and
y values. This gives the expression
χ2pos =

1
N

N

˜ i ) − x i ]2 + [ y(z
˜ i ) − y i ]2
[x(z
σ2pos

i =0

(3.29)

where N is the number of data points in the experimental track, and σpos is an expected deviation. We
chose σpos to be 0.5 cm based on the average size of the pads.
The energy component of the objective function is produced by projecting each generated track onto
the pad plane and comparing the resulting hit pattern to the experimental hit pattern. This gives a function of the form
χ2en =

∆A 2
1
.
Nhit hit pads σ2en

(3.30)

Here A represents the total charge deposited in a particular pad. The estimated uncertainty σen was
chosen to be 10 % of the maximum value of A present in the experimental track. The sum is taken over
all pads that were hit by the experimental track, and Nhit is the number of such pads. Limiting the sum
to pads present in the experimentally measured event helps the fit converge when parts of the track are
missing since otherwise χ2en would have an unbalanced contribution from all the pads that are present
in the generated track but not the experimental one.
45

The final component of the objective function helps constrain the reaction vertex to be near the z
axis. Normally, we could fit the beam track and ensure that the beam track and the track of the scattered
particle intersect, but in the 46 Ar data set, the beam tracks were not usable (see Section 5.3.5 for details).
This can be addressed by adding the following extra term to the objective function:
χ2vert =

x 02 + y 02
σ2vert

.

(3.31)

This is a simple quadratic constraint on the (x, y) location of the vertex that imposes a penalty on the fit if
the vertex is not near the axis. One problem with this approach is that it assumes that the beam is exactly
along the z axis, which it might not be due to emittance; however, the emittance of the beam seems to
be small enough that this approximation is acceptable.
The complete objective function is, then, the sum of Eqs. (3.29), (3.30), and (3.31):
χ2 = χ2pos + χ2en + χ2vert .

(3.32)

3.4.4 Finding initial parameter values for optimization
The optimization algorithm must be given a starting point before it can begin generating candidate
tracks. The starting point is found using a simple linear fit. In Section 3.3.3, we found that the arc length
coordinate s = ρθ can be plotted as a function of the z coordinate to get a straight line. This was also
shown in Fig. 3.7a. This technique can also be used to get a seed for the Monte Carlo algorithm.
To find the seed, we can use the previously found center of curvature to find ρ and φ using the equations in Eq. (3.24). This time, however, we subtract the angular coordinate of the first point in the track to
make φ start at 0. ρ and φ are then multiplied together to get the arc length coordinate s, and a line is fit
to the plot of s vs z. The slope of this line defines the initial polar angle θ0 of the track, and the intercept
of the line defines the z 0 position of the reaction vertex.
The initial energy of the particle is calculated from the radius of curvature. In a magnetic field B , a
charged particle with transverse velocity v, mass m and charge q will move in a circle of radius
ρ=

mv
.
Bq

(3.33)

This particle’s kinetic energy is, of course,
1
E = mv 2 .
2
46

(3.34)

Combining these yields
E=

B 2ρ2 q 2 B 2ρ2 Z 2e 2
=
,
2m
2Am p

(3.35)

where A and Z are the mass and charge number of the particle. This can be used to find an estimate for
the initial energy E 0 , assuming we know the identity of the particle.
Finally, the remaining coordinates x 0 and y 0 are assumed to be zero since the data calibration procedure places the beam along the z axis.

3.5 The Monte Carlo optimizer
3.5.1 Description of algorithm
Tracks from the AT-TPC can be fit using a simple Monte Carlo algorithm. The algorithm starts with the
seed point found above. Then, it generates a large number T of candidate parameter sets from a uniform
distribution over the six-dimensional (x 0 , y 0 , z 0 , E 0 , φ0 , θ0 ) parameter space. The distribution is centered
on the seed point in each dimension, and the width is set to a configurable value. Each of these candidate
tracks is then simulated, and the values of χ2 are calculated using the equations in Section 3.4.3. The
track with the lowest value of χ2 is selected, and the parameter space is re-centered around this new
point before being compressed in each dimension by a multiplicative reduction factor R. This process is
then repeated for a fixed number of iterations I . At the end of the I iterations, the best track is accepted
as the fit result. Figure 3.10 shows how the parameter space converges to a solution across the iterations
of this algorithm.
From the description of this algorithm, it should be apparent that this is not the most efficient way
to fit a track. For each event in the data, the computer has to essentially guess the answer N = I T times
and pick the best result. Fortunately, the convergence of the size of the parameter space helps the algorithm move toward a solution, and it prevents it from generating candidate tracks that are too far from
the minimum. Furthermore, although generating such a large number of tracks is certainly computationally costly, modern high-performance computers and parallel programming techniques help reduce
the runtime to something reasonable.

47

Figure 3.10: Samples drawn by the Monte Carlo algorithm while fitting one event. The horizontal axes
represent the sample number, while the vertical axes represent the value of each of the 6 variables. x, y,
and z are given in meters, E is in MeV, and θ and φ are in radians. Each 500-sample iteration of the code
is visible as a large-scale block in this plot.
3.5.2 Checking convergence
The Monte Carlo algorithm described above does not include an explicit test for convergence, so the
parameters of the fitter must be chosen in such a way that convergence is likely. There are three main
parameters that control the convergence of the fit: the number of candidate tracks T generated per iteration, the number of iterations I , and the reduction factor R by which the parameter space is compressed
after each iteration.
The values of T and R must be chosen carefully to prevent the fitter from converging to local minima.
Much like in the process of simulated annealing, it is important to reduce the size of the parameter space
very slowly and to sample as many points in the space as possible on each iteration. The main limitation
on T is the available computing power, and the value of R must be balanced with the value of I to control
the width of the parameter space at the end of the fit.
If the parameter space has a width ∆X 0 along dimension X at the beginning of the fit, the final width
in that dimension will be
∆X = (∆X 0 )R I .
48

(3.36)

Reduction factor

0.9 0.8 0.7 0.6 0.5

100 points
13
18
29
27
1

17
17
24
41
17

17
21
27
47
31

16
25
30
45
60

200 points
15
19
39
51
64

10 15 20 25 30

24
34
38
45
8

25
37
54
56
26

26
26
42
61
59

22
29
47
64
71

300 points
27
27
46
69
77

10 15 20 25 30

22
40
47
52
6

22
43
54
67
34

31
41
63
72
68

32
36
57
77
75

400 points
26
42
62
76
79

10 15 20 25 30

38
43
51
54
10

23
53
64
71
48

34
46
62
69
72

40
41
65
83
79

500 points
28
37
60
84
84

10 15 20 25 30

39
48
58
57
14

40
56
68
76
46

34
53
69
78
74

33
44
66
87
78

34
44
76
90
86

10 15 20 25 30

Number of iterations Number of iterations Number of iterations Number of iterations Number of iterations
Figure 3.11: Results of the Monte Carlo convergence study. Subplots represent different numbers of
points T per iteration. Within each plot, rows represent the reduction factor R, and columns represent
the number of iterations I . The number in each cell represents the percentage of fits in which the optimizer converged to the correct minimum, subject to a tolerance of 10 mm in each dimension for the
vertex position, 100 keV for the initial proton energy, 5° for the azimuthal angle, and 2° for the scattering
angle. The fit was run 1000 times on the same simulated track for each set of parameters.
I

T

R

σx,y,z

σE

σφ

σθ

20

500

0.8

0.1 m

4.0 MeV

60°

30°

Table 3.1: Monte Carlo fit parameters I , T , and R and initial parameter space widths σ used in the
46
Ar(p, p ) measurement.
This gives an estimate for the uncertainty of the fit result in each dimension. This final uncertainty can
be reduced by either using a larger value of I or a smaller value of R.
The values of these parameters can be optimized by fitting one event repeatedly and examining the
distribution of the fit results. The fit parameters can then be varied to find the best choices. This was
done for the 46 Ar(p, p ) data set by simulating one proton track and then fitting it 1000 times for various
combinations of parameters. The results of this study are shown in Fig. 3.11, and the parameters used in
the analysis are shown in Table 3.1.

3.6 Alternative fitting algorithms
Several other minimization algorithms were considered in addition to the Monte Carlo optimizer described above. Although each of these techniques has some advantages, the naïve Monte Carlo algorithm
ultimately produced the best, most consistent results. For completeness, a few of the other algorithms
we tried are described in this section.

49

Approx. first derivative

Objective function

400
300
200
100
0
100
200
300
400
500
8
7
6
5
4
3
2
1
0

1.6

1.8

2.0

2.2

2.4

Energy [MeV]
(a) Surface

(b) Slice of surface

Figure 3.12: Simulated two-dimensional error surface and a one-dimensional slice through that surface.
These were generated from a simulated 46 Ar(p, p ) scattering event by varying the scattered proton’s initial energy and scattering angle, calculating the objective function at each point. 400 points were simulated in each dimension. On the surface plot, the horizontal axis represents energy in MeV, the vertical
axis represents the energy component of the objective function, and the axis going into the page represents the scattering angle in radians.
3.6.1 Gradient-based methods
The default first approach to an optimization problem is often to use a gradient-based technique. These
algorithms are generally robust and computationally efficient [28], and commonly used software libraries like SciPy [25] and

MINUIT

[24] implement a variety of these algorithms for easy use. However,

these techniques are not well-suited to the problem of fitting tracks in the AT-TPC since smooth derivatives are not available for our objective function.
As described previously, the value of the objective function for a given set of track parameters is found
by simulating a track in the AT-TPC and projecting it onto the pad plane. This produces an objective
function that is not smooth when examined closely, which causes derivatives calculated using the finite
differences method to diverge.
Kolda et al. [28] state that simulation-based optimization problems like ours often produce objective functions that are nonsmooth and sometimes discontinuous. They name conditional branching
(if/then/else constructs), convergence tests, adaptive algorithms, and loss of numerical precision as a
few common causes of nonsmoothness in simulation-based problems. They also point out that the com-

50

Algorithm 3 Basic compass search algorithm [28]
Given a function f : Rn → R, an initial guess x 0 ∈ Rn , a tolerance ∆tol for determining convergence, and
an initial value ∆0 > ∆tol for the step-length control parameter,
function C OMPASS( f , x 0 , ∆tol , ∆0 )
for k = 1, 2, 3, . . . do
D ⊕ ← {±e i |i = 1, 2, . . . , n}, the set of unit vectors spanning Rn
if there exists a d k ∈ D ⊕ such that f (x k + ∆k d k ) < f (x k ) then
x k+1 ← x k + ∆k d k
Take step in best direction
∆k+1 ← ∆k
No change to step size
else
x k+1 ← x k
No step is taken
1
Shrink step size
∆k+1 ← 2 ∆k
if ∆k+1 < ∆tol then
return x k+1
end if
end if
end for
end function
plexity of the simulation code often makes it difficult or impossible to find an analytic representation of
the derivative, even when automatic differentiation libraries are available.
Figure 3.12 shows the energy component of the objective function evaluated for a simulated event as
a function of energy and scattering angle. The surface shown in Fig. 3.12a is quite rough when examined
on a fine scale. This is more apparent in Fig. 3.12b, which shows a slice taken through the error surface at
the scattering angle that corresponds to the minimum. The derivative of this slice with respect to energy,
as estimated by finite differences, is very rough and contains abrupt jumps. This would cause problems
for any gradient-based optimization technique.

3.6.2 Compass search
Kolda et al. [28] describe an early, basic optimization algorithm for nonsmooth functions that they call
compass search. This technique, which is described fully in Algorithm 3, evaluates the objective function
at a given initial point and then searches along a set of direction vectors for a better point. If a better
point is found, the fitter moves to that point. If no better point is found, the fitter shrinks its step size to
prevent it from stepping over the minimum.
One problem with this algorithm is that it is very slow to converge if the direction from the initial
point to the minimum is not along one of the unit vectors e i . In this case, the fitter will “bounce” between

51

OPTIMIZATION BY DIRECT SEARCH

429

Figure 3.13: An illustration of the points tested in the Nelder-Mead simplex method. The vertex v 2 is
(a) across the point c using four different values
(b) of the multiplicative factor α, leading to the four
reflected
possible new vertices, shown in purple, yellow, blue, and green. (Reproduced from [28, Fig. 4.2(b)].)

. 4.2 Examples in R2 of both (a) the single possible step for the simplex algorithm of Spendley,
Hext, and Himsworth
and (b) theinfour
basic pattern,
steps forevaluating
the Nelder-Mead
simplex
algorithm.
two of the directions
a zig-zag
the objective
function,
which might be computa-

tionally expensive, many more times than it needs to. This could be corrected by allowing the fitter to

try a linear combination
of direction
vectors
add this new
direction
The simplex algorithm
of Neider and
Mead [194]
is and
a variation
on this
basicvector
idea to the set of e i for later
at allows, in effect, a simple line search of the form vn + a(c ? vn), with a set
iterations, but this causes the set of e i to contain more vectors than are needed to span the parameter

four possible choices for a. Typical choices are a G {^,?,2,3}, as illustrated
n
Figure 4.2(b);space
see R
Lagarias
et al.
[164]
for a particularly careful and complete
, which also
leads
to inefficiency.

scription of the Nelder-Mead simplex algorithm. The line search has the effect of
owing the shape of the simplex to deform (for any choice of a other than 2), which
touted as a feature
allows the
simplex
to adapt to the local topology of the
3.6.3 that
Nelder-Mead
simplex
algorithm
nction, hence references to this algorithm as the adaptive simplex method.
While these two
simplex
algorithms
arethat
classical
direct
methods,
neither above
is
A more
sophisticated
option
is similar
to thesearch
compass
search described
is a simplex method,
GSS method. These two simplex methods search along the single search direction
such as the algorithm described by Nelder and Mead [39]. In this method, rather than evaluating the
? vn). Further, both these methods enforce only simple decrease, but in a sense
at is subtly different
from the
definition
of simple
decrease
defined
inin(2.3)
and
objective function
along
a set of directions
e i from
the initial
point, as
the compass
search, the function
d in Step 2 of Algorithm 3.2, since their step acceptance condition
requires
simple
is evaluated at the vertices of a simplex of n + 1 points in Rn . For example, in two-dimensional space R2 ,
rease in / at the vertex in the simplex with the second highest function value vn_i,
t at the vertexthe
infunction
the simplex
with at
the
function
value
(denoted
is evaluated
thelowest
verticesknown
of a triangle
(a set of
n+1 =
3 points).vo).
On each iteration, the vertex
us both algorithms only ensure improvement of the function value at the sequence
with the largest value of the objective function is reflected across the center of the opposite face of the
worst vertices,
but it is the sequence of best vertices (i.e., those with the lowest
nction value) simplex,
that ultimately
of interest. As an interesting aside, the Nelder
c = n1 (v 0 +vis
1 +· · ·+v n ) [28]. The function is then evaluated at a series of points v r = v i +α(c −v n )
ad simplex algorithm in R1 can be restated as a GSS method. Since in this special
3
for as
a predetermined
values
α, commonly
α ∈ 12 ,with
[28].next
The new vertex v r with the
2 , 2, 3 the
e vo serves both
the centroid set
of of
the
opposite
face andchosen
as theasvertex
ghest function smallest
value, in
effect
searchfunction
is alongis two
(a value
positive
basis
value
of thethe
objective
then directions
accepted if its
is less
thaninthe value of the function
from ^o, and any improvement in the function value is with respect to / at vq.
at thefor
vertex
the second-highest
value of the
objective function
(The second-worst vertex is
alternative proof
thiswith
special
case, under different
assumptions
on /, [28].
is given
[164].
chosen as the point of comparison since it will then be the worst vertex in the next step [28].)
McKinnon [179] constructed a family of functions in R2 which demonstrates that
This algorithm
tends
to fail
perform
better thantothe
search
method
e Nelder-Mead simplex
algorithm
can
to converge
a compass
stationary
point
of /,since it evaluates the oben if the family
is parameterized so that / is strictly convex and has up to three
jective function at fewer points on each iteration; however, McKinnon [35] showed that there exists a
ntinuous derivatives. Key to the failure McKinnon demonstrates is the ability to
familyRepeated
of degenerate
functions that
cause
thethe
simplex
to collapse
along the direction
of steepest descent,
form the simplex.
deformations
can
cause
sequence
of simplices
pro
ced by the Nelder-Mead simplex algorithm to converge to a degenerate simplex.
preventing further progress. McKinnon refers to this phenomenon as “repeated focused inside contrac-

tice in Figure 4.2 that choosing a = \ (a so-called inside contraction) replaces
with a vertex that moves closer to c. In McKinnon's examples, the Nelder-Mead

mplex algorithm repeatedly chooses a = \ (only) and the52simplices converge to a
raight line that is orthogonal to the steepest descent direction and have interior
gles which tend to zero (i.e., the simplices collapse along the steepest descent di

152

K. I. M. MCKINNON
(0)

v

Downloaded 04/08/16 to 35.9.62.246. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

y

(2)

v

x

(3)

v

(1)

v

Fig. 2.1. Successive simplices with rfics.

Figure 3.14: Steps taken by a vertex during simplex collapse, a failure mode of the Nelder-Mead algorithm. (Reproduced from [35, Fig. 2.1].)
shown in Figure 2.1. The general form of the three points needed at one step of the
Nelder–Mead method is therefore
tion,” or RFIC. This is illustrated in Fig.
3.14. Although our objective function does not necessarily cause
(2.7)
v (n) = ( n1 , n2 ),
n+1
this behavior,
one way
(2.8)it serves to illustratev (n+1)
= ( in
, n+1
), method can fail.
1 which
2 this
(2.9)

r(n) = (

n
1 (1

1 ),

n
2 (1

2 )).

Provided
(2.1) and (2.3) hold at these points, the simplex method will take the
3.6.4 Simulated
annealing
inside contraction step assumed above.
Note that the x coordinates of v (n) and v (n+1) are positive and the x coordinate
(n)
Unlike the of
deterministic
methods described previously in this section, simulated annealing is a stochasis negative.
r

ˇ
. Consider
f (x,ety)al.
given
3. process.
Functions
which
cause RFIC
tic optimization
It was
independently
described
[5]the
by function
Kirkpatrick
[26] by
in 1983 and Cerný
2
(3.1)
, xthe
 0,
f (x, y) by
= ✓comparison
|x|⌧ + y + ywith
[10] in 1985.
This algorithm was created
physical process of annealing in sta⌧
2
= ✓x + y + y , x 0,
tistical mechanics. The basic procedure consists of creating a Markov chain of candidate solutions by
where ✓ and are positive constants. Note that (0,-1) is a descent direction from
drawing from
a random
accepting
new states
a probability
the origin
(0,0)distribution
and that f isand
strictly
convexthe
provided
⌧ > with
1. f has
continuousthat
firstdepends on
derivatives if ⌧ > 1, continuous second derivatives if ⌧ > 2, and continuous third
a “temperature” parameter, which is slowly cooled. A basic description of the algorithm, adapted from
derivatives if ⌧ > 3. Figure 2.2 shows the contour plot of this function and the first
of the
method for
Bertsimas two
and steps
Tsitsiklis
[5],Nelder–Mead
is given in Algorithm
4. the case ⌧ = 2, ✓ = 6, and = 60. Both
steps are inside contractions.
One thing that is immediately apparent about this algorithm is its similarity to the simple Monte

Carlo algorithm described in Section 3.5.1. In each case, we look for a new point on each iteration that

53

Algorithm 4 Basic simulated annealing
Given a function f : Rn → R, an initial state x 0 ∈ Rn , and a decreasing set of temperatures Ti
function A NNEAL( f , x 0 , Ti )
for i=1, 2, . . . , N do
x ← x i −1 + a random step in some direction
if f (x ) ≤ f (x i −1 ) then
Accept x as the new value x i
else
Accept x as the new value x i with probability e −(x −x)/Ti .
Otherwise, x i ← x i −1 .
end if
end for
return x N
end function
satisfies a criterion that we slowly make more stringent. The only concrete differences are that in simulated annealing, the next point is chosen by taking a random step from the current point, and sometimes
we accept a new point with a higher value of the function being minimized.
Computationally, this algorithm can be more efficient than the naïve Monte Carlo method since it
will ideally test fewer points on each iteration. However, if the temperature parameter is too small, the
algorithm might not be able to quickly find a better point than the current one, leading it to test a large
number of points on each iteration. This can quickly make simulated annealing slower than the naïve
Monte Carlo algorithm.
Ultimately, simulated annealing provided no clear benefit over the naïve algorithm for data from the
AT-TPC. It did not converge reliably, and the resolution was poorer. Finally, although simulated annealing
is a more efficient algorithm in principle, the predictable run time of the simple Monte Carlo method is
attractive when scheduling the analysis on a high-performance computer.

3.6.5 Neural networks
Modeled after networks of neurons in the brain, an artificial neural network consists of a set of interconnected nodes where each connection is assigned a weight [8]. The weights are chosen such that when a
set of inputs is applied on one side of the network, the correct outputs are produced on the other side. The
network is “trained” with a large data set with known inputs and outputs by presenting the inputs one by
one and adjusting the weights on each iteration. The goal is to generalize the network’s “understanding”
of the problem so that it can successfully find the correct outputs for an unfamiliar set of inputs.

54

Neural networks were considered as a potential first step to the fitting process for AT-TPC data. However, the technique using the Hough transform and linear fits described in Section 3.4.4 does a much
better job and is much simpler.

3.6.6 Kalman filter
The Kalman filter is a recursive optimization algorithm that is commonly used in diverse fields such as
automated navigation, robotics, and aviation [30, 56]. It also finds applications in high-energy physics
[14], making it a good candidate for fitting tracks in the AT-TPC.
A Kalman filter models a system using a state vector and a physical model. The state vector represents
the state of the system at a given moment, and might include components like the position of a particle
and its momentum. The physical model must be able to take the state vector and produce an estimate
of the new state vector at the next moment in time. For the AT-TPC, a suitable model might include the
Lorentz force and a contribution to account for the stopping power of the gas.
Rather than considering the data set as a whole, the Kalman filter fits each point individually, iterating
through the data set only once. On each step, the model is used to predict the new state vector based on
the state from the previous step. The measured data point is then averaged with this prediction to correct
it. This averaging allows the Kalman filter to smooth out noisy data [30].
One drawback of most implementations of the Kalman filter is that they require data points that are
evenly spaced in time. This is very often not the case in the AT-TPC, however, since there are frequently
large gaps in the tracks recorded by the detector. This limitation ultimately prevented us from using the
Kalman filter to obtain precise results.

3.7 Vertex energy reconstruction
Using the fit results, we have enough information to reconstruct the energy of the beam projectile at
the reaction vertex in two different ways. This redundancy is a benefit since it provides a way to check
whether the analysis is correct, and it also provides a measure of the resolution of the detector. Furthermore, any systematic deviation between the two reconstructions can be used to separate elastic and
inelastic scattering events since the position of the vertex along the z axis does not depend on the kinematics.

55

3.7.1 Using the vertex position
The Monte Carlo result gives the location of the vertex within the detector. The energy of the beam particle at the vertex can be calculated by finding the amount of energy lost by the beam particle to the gas
as it moved from the beam entrance window to the vertex.
Assuming the beam particle entered the chamber with energy E i , the stopping power calculation
from SRIM [57] gives the particle’s range as a function of energy, R(E ). Thus, the initial range of the fullenergy beam particle is
R i = R(E i ).

(3.37)

The beam particle loses energy as it travels from the window to the vertex. Since R is a function of energy,
we could also think of this as the particle losing range. Then, the amount of range “remaining” at the
vertex is
R f = R i − (1.0 m − z 0 ).

(3.38)

Finally, the inverse of the range function, R −1 , can be used to find the energy of a particle with range R f :
E f = R −1 (R f ).

(3.39)

This final energy E f is the energy of the beam particle at the reaction vertex.

3.7.2 Using kinematics
After fitting the tracks, the initial energy and scattering angle of the recoil particle are known. Assuming that this particle was at rest before it interacted with the projectile, we have enough information to
reconstruct the vertex energy using kinematics.
In the laboratory frame, the target nucleus has four-momentum P tlab = (m t , 0, 0, 0). In the center-ofmass frame, its four-momentum is P tCM = (E t , −p CM , 0, 0). These four-vectors can be used to find the
components of the Lorentz boost matrix from the center-of-mass frame to the lab frame. The results are

cosh χ =

2
p CM
+ m t2

mt
p CM
sinh χ =
,
mt

56

(3.40)
(3.41)

where χ ≡ tanh−1 (v/c) is the rapidity. The Lorentz transform can then be applied to the recoil particle to
find p CM . The lab energy of the recoil particle is
E rlab = E rCM cosh χ − p CM cos θCM sinh χ.

(3.42)

2
If we plug in Eqs. (3.40) and (3.41), assume elastic scattering, and solve for p CM
, we find
2
p CM
=

m t E rlab − m t2
1 − cos θCM

.

(3.43)

Finally, the kinetic energy of the beam particle Tb can be found using the relativistic invariant
2

s=

(3.44)

Pi
i

for four-momenta P i . We can find this invariant for the two frames, set them equal to each other, and
solve for Tb to find
Tb =

1
2m t

2
m b2 + p CM
+

2
m t2 + p CM

2

− (m t + m b )2 ,

where m t is the mass of the target and m b is the mass of the beam.

57

(3.45)

CHAPTER 4
SIMULATIONS

In order to understand the behavior of the AT-TPC and establish expectations for experimental results,
we simulated several sets of data under the same conditions as the experiment. The purpose of this was
to model the resolution of the detector as well as the efficiency of the trigger and analysis methods. By
modifying the number of simulated events, we can also determine how many events are needed in an
experiment to be able to resolve resonances and measure their strengths.
The simulations discussed in this chapter were performed using software that is based on the same
libraries as those used for analyzing AT-TPC data. Specifically, the simulated particle tracks were generated using the same functions that are used to generate candidate tracks in the Monte Carlo algorithm,
so the simulated tracks use the same energy loss data, Lorentz force calculations, calibrations, and pad
plane projection functions as the analysis process. Therefore any differences between the simulations
and the data should only be due to features of the data or effects neglected in the model, and not due to
discrepancies between the simulation and analysis codes.
This chapter begins with an overview of the simulation process. Following that are the results of
several simulations that were performed.

4.1 Overview of the simulation process
Although several simulations with slightly different parameters and settings were performed, the simulation process was largely the same for each of them. This process is described in general terms in this
section; any deviations from this procedure for specific simulations are outlined in later sections.

4.1.1 Parameter generation
The first step of the simulation is to create a set of simulated track parameters. This can be done in a
number of different ways depending on the goals of a particular simulation.

58

v0
y
b

θ0
y0

window
L
θB

z0

focus

z

Figure 4.1: Illustration of coordinate calculation for uniform distributions. The meanings of the variables are given in Section 4.1.1.1. Note that this is a projection into the z y plane, and the true threedimensional orientation also depends on azimuthal angles φ0 and φB and vertex location x 0 .
4.1.1.1 Uniform distribution
The first set of simulations was performed using parameters drawn from flat distributions. This is useful
for modeling the efficiency of the detector since it uniformly populates the parameter space of the tracks,
making it easy to see what angles and energies are cut by the trigger or by the analysis. Flat distributions
were generated as follows:
1. Pick z 0 from a uniform distribution between 0 mm and 1000 mm, which are the locations of the
sensor plane and the beam entrance window, respectively.
2. Generate a random beam vector b by selecting the beam azimuthal angle φB and beam polar angle θB such that the result will be uniformly distributed over a portion of the unit sphere. The
azimuthal angle was allowed to vary uniformly over [0, 2π), and the polar angle varied between
0° and 2°, where 0° is parallel to the beam axis. This beam vector helps simulate the effect of the
emittance of the beam.
3. Use the beam vector to calculate the values of x 0 and y 0 at the vertex location z 0 . Assume that the
beam tracks originate at a focus 14 cm upstream of the beam entrance window. The location of
this focus was chosen such that a particle leaving it at the maximum beam polar angle of 2° would
make it through the beam entrance window, which had a diameter of 10 mm.
4. Generate the azimuthal and polar angles describing the outgoing proton’s initial velocity vector v0
with respect to the beam track. Select the azimuthal angle from a uniform distribution over [0, 2π]

59

and the polar angle from a uniform distribution over [π/2, π). The polar angle is greater than π/2
since the beam is moving in the −zˆ direction.
5. Transform the two angles from the previous step to find angles with respect to the coordinate axes
rather than the beam track. This transformation corresponds to a rigid-body rotation, and it is
done with Euler angles. This yields the proton track’s initial azimuthal angle φ0 and polar angle θ0 .
6. Find the length L of the beam track between the entrance window and the reaction vertex. Use
this to calculate the 46 Ar vertex energy with the procedure from Section 3.7.1. Calculate the proton
energy E 0 using the kinematics formulas below. If the proton energy is less than 0.5 MeV, discard
the event as it will almost certainly be fit poorly.
The relationships between the different quantities described above are illustrated in Fig. 4.1.
The initial energy of the proton is calculated from the 46 Ar vertex energy using kinematics formulas.
To derive these equations, first consider a generic scattering process B (A, D)C where beam particle A
has kinetic energy T and target B is at rest. Let C be the target-like recoil particle and D be the beamlike ejectile. Since we are interested in finding the proton energy, we will need to calculate the energy of
particle C . In the center-of-mass frame, this is
ECCM =

p CM

2

+ mC2 ,

(4.1)

where p CM is the outgoing momentum in the center-of-mass frame. This momentum is
p CM =

2
s − mC2 − m D

4s

2

2
− 4mC2 m D

(4.2)

where s is the square of the total incoming four-momentum
s = (P A + P B )2 = (m A + m B )2 + 2m B T,

(4.3)

which is a relativistic invariant.
Similarly, in the lab frame, the energy of particle C is
EC =

pC2 + mC2 .

(4.4)

The lab momentum pC can be found be equating the components of the lab and center-of-mass momenta that are transverse to the Lorentz boost between the two frames:
pC =

p CM sin θCM
sin θC
60

.

(4.5)

Since the lab scattering angle θC is known, this leaves the center-of-mass angle θCM to be determined.
This can also be found by equating components of the two momenta, yielding
sin θCM =

sin θC ECCM cos θC sinh χ + α

p CM cosh2 χ − cos2 θC sinh2 χ

(4.6)

where
α≡

cosh2 χ ECCM

2

+ mC2 − cosh2 χ + cos2 θC sinh2 χ .

In these equations, χ is the rapidity of the boost from the center-of-mass frame to the lab, which is given
by


χ = log 

p CM +

2
m B2 + p CM

mB



.

(4.7)

The initial center-of-mass frame momentum p CM can be found by using Eq. (4.2) and replacing mC and
m D with m A and m B , respectively. These equations are evaluated with the appropriate mass values to
find the initial energy of the proton in the lab frame from the randomly generated scattering angle and
46

Ar vertex energy.

4.1.1.2 Simulated scattering distribution
Although the uniform distributions described above are useful for modeling the efficiency of the detector, they are not representative of the distributions of parameters encountered in a real experiment. A
more realistic set of simulated tracks can be created by drawing parameters from a distribution created
by a simulation of the scattering process to be measured.
For the simulation of the 46 Ar(p, p ) experiment, the parameter generation was based on a calculation
from the R matrix code

DSIGMAIV 1

using the expected resonance properties shown in Table 4.1. The

cross section of the reaction was calculated using this code for center-of-mass incoming proton energies
from 2.0 MeV to 4.5 MeV and center-of-mass scattering angles from 40° to 140°. The elastic scattering
component was simulated using the global optical potential from Koning and Delaroche [29]. The optical
potential parameters used are shown in Table 4.2 and the results of the calculation are shown in Fig. 4.2.
The results from DSIGMAIV give the differential cross section as a function of proton energy and scattering angle in the center of mass frame. The cross section is effectively an unnormalized probability of
1 Originally written by Wang Hongwei of the Shanghai Institute of Applied Physics at the Chinese Academy of Sciences. The
version used included modifications made by Wolfgang Mittig of the NSCL.

61

47

Ar state
g.s.
1.130

ER

L

2J

Π

S

Γp

2.582
3.712

1
1

3
1

−1
−1

0.6
0.8

10.94
63.80

Table 4.1: Expected resonance properties calculated from the results of Gaudefroy et al. [18] using the
procedure from Section 1.1.8. The resonance energy E R is given in MeV, and the resonance width Γp ,
as calculated by DSIGMAIV, is given in keV. L is the orbital angular momentum, J is the total angular
momentum, Π is the parity, and S is the spectroscopic factor.

Term

Vr

Rr

ar

Vi

ri

ai

Volume
Surface
Spin-Orbit

62.60
–
5.80

1.19
–
1.00

0.67
–
0.59

0.30
7.70
–

1.19
1.29
–

0.67
0.54
–

Table 4.2: Global optical model [29] parameters used in the simulation. Vr , R r , and a r are the depth,
radius, and diffuseness parameters for a real-valued Woods-Saxon potential, while Vi , R i , and a i are the
corresponding parameters for an imaginary Woods-Saxon potential. The depth is given in MeV, and the
radius and diffuseness are given in fm. The Coulomb radius used was 1.27 fm.

120
103

Scattering angle [deg, CM]

Differential cross section [mb/sr]

140

100
80

102

60
2.0

2.5

3.0

3.5

Proton energy [MeV, CM]

4.0

4.5

Figure 4.2: Results of DSIGMAIV calculation in the region of interest.

62

observing a particle with the given energy and scattering angle, so these results can be used with rejection sampling to generate a set of proton energies and scattering angles that follow our expectations for
the experimental results.
Given the incoming proton energy E pCM and center-of-mass scattering angle θCM , we can calculate
the energy of the outgoing proton using relativistic kinematics. The Lorentz boost matrix from the centerof-mass frame to the lab frame is given by

cosh χ sinh χ


 sinh χ cosh χ
Λ=

 0
0


0
0


0 0


0 0


1 0


0 1

(4.8)

for rapidity χ. The values of sinh χ and cosh χ can be found by boosting the target proton from the centerof-mass frame to the lab frame, where it is stationary. The result of that is
2
p CM
+ m p2

cosh χ =

(4.9)

mp

sinh χ =

p CM
,
mp

(4.10)

where m p is the proton mass and p CM is the momentum in the center-of-mass frame. The latter quantity
is known since we know E pCM :
E pCM

p CM =

2

− m p2 .

(4.11)

Finally, we can boost the outgoing proton’s four-momentum from the center-of-mass frame to the lab
frame to find E plab :
E plab = E pCM cosh χ − p CM cos θCM sinh χ
=
=

2
p CM
+ m p2
2
p CM

mp

2
p CM
+ m p2

mp

− p CM cos θCM

(1 − cos θCM ) + m p .

p CM
mp
(4.12)

After subtracting away the rest mass, this value can be used as the initial energy E 0 in the parameter set.
The scattering angle in the lab frame is given by
θplab ≈

1
π − θpCM .
2
63

(4.13)

This is a non-relativistic approximation, but it holds quite well at the low energies of this experiment.
This value is used to find the parameter θ0 , which is given by
θ0 = π − θplab

(4.14)

since the coordinate system is defined such that the beam moves in the −zˆ direction. The azimuthal angle φ0 is a free parameter in the scattering process, so it is chosen randomly from a uniform distribution
over the domain [0, 2π).
Next, we need to calculate the initial position of the particle. The transverse positions x 0 and y 0
are chosen randomly from normal distributions centered at 0 with standard deviation 10 mm, a choice
which simulates the emittance of the beam. The z 0 position is found using the kinematics of the beam
particle and the stopping power of the gas.
The energy of the 46 Ar beam particle in the center-of-mass frame is simply
CM
E Ar
=

2
2
p CM
+ m Ar

(4.15)

for the same p CM as in Eq. (4.11) above. The lab energy is then
lab
CM
E Ar
= E Ar
cosh χ + p CM sinh χ.

(4.16)

Expressions for p CM , cosh χ and sinh χ from above can then be substituted into this equation to find the
value.
Once the vertex energy of the 46 Ar nucleus is known, the vertex position z 0 can be found from the
calculated range of the particle in the detector fill gas, a quantity that can be calculated with SRIM [57].
If the range of the incoming 46 Ar nucleus at the full beam energy is R i and the range of an 46 Ar nucleus
with the final vertex energy is R f , then the distance traveled by the particle must be R i −R f , so the vertex
is located at
z 0 = (1.0 m) − R i − R f .

(4.17)

At this point, we have a complete set of parameters (x 0 , y 0 , z 0 , E 0 , φ0 , θ0 ) that describes the simulated
proton tracks.

4.1.2 Tracking and event generation
Once the set of initial parameters has been found, the proton tracks can be simulated. The simulation
method is the same as in Section 3.4, but instead of finding the peak in each channel, each entire signal
64

is kept. These signals represent the ideal response of the detector. To better simulate the behavior of the
detector, three kinds of noise are added to the simulated signals.
First, random Gaussian noise is added to each signal to simulate noise in the electronics. The distribution used to generate this noise has a mean of 0 and a standard deviation of 6, where both values
are measured in ADC bins. This standard deviation was found by analyzing the distribution of the noise
amplitudes of a number of experimental signals inside a time window where the traces contained no
pulses.
The second kind of noise added is a simulated baseline depression signal of the form discussed in
Section 3.1. The track of an 46 Ar nucleus with the full beam energy is simulated in the detector, and the
signals are accumulated bin-wise to find the total energy deposited as a function of drift time. This signal
is then inverted, rescaled, and added to each simulated signal to approximate the baseline depression
effect in the experiment. The rescaling amplitude was found empirically from recorded experimental
signals.
Finally, a constant offset, or pedestal, is added to each channel. The pedestals are slightly different for
each channel as they originate from the baseline of each channel’s ADC. These pedestals were measured
during the experiment by triggering the detector with a pulser. After adding the pedestals, the signals are
clipped to ensure they remain between 0 and 4096, the range of the 12-bit ADC.

4.1.3 Simulating the trigger
One of the most complex parts of the simulation process is determining whether an event would trigger the GET electronics. Each channel in the AGET chip has a subsystem that compares the shaped,
amplified input signal to a configurable threshold to determine if the channel was hit (Fig. 4.3). These
channel-level triggers are then accumulated at the CoBo level to make a decision about whether an event
should be recorded. This process must be simulated in order to estimate the trigger efficiency during the
experiment.
All of the information about the GET electronics presented in this section was taken from Baron and
Delagnes [3].

65

The output signal of the common filter (Fig. 8) is amplified through a
differential buffer to cover an input dynamic range of the discriminator. This range
can be fixed to 5% or 17.5% of input dynamic range of the analog channel. The
differential gain is controlled by the bit state2<24>.

SCA
CSA

PZC

FILTER

G2
512 cells
Trigger
+
- Discri

Gain

Sign

Polarity
3-bit
DAC

4-bit
DAC

Polarity

Common to all channels

Hit
channel
register

4

inhibit

Slow Control Register
[channel level]

Slow Control Register
[Asic level]

Fig. 8: block schematic of the channel.

Figure 4.3: Schematic of a channel in the GET electronics, showing the trigger threshold circuit. (Reproduced fromThe
[3].) threshold value of the discriminator is controlled and adjusted through 2

DACs. The first is global to all the channels. Its 3 bits give the MSB of the
threshold value plus 1 bit for polarity (state1<22:19>). The second is attached to
4.1.3.1 Generating channel-level trigger signals
TheGET-QA-000-0005
trigger threshold for each channel can be represented as a 7 bit integer. The three most-significant
AGET: Data Sheet. Production Version 2.0. Document 1.0

13/90

bits are common to all channels in a given AGET, and the four least-significant bits are unique to each
channel. An additional bit is used to store the polarity.
This raw, integer threshold value must be converted to be on the same scale as the data. The trigger
discriminator in the electronics covers 17.5 % of the input range of the ADC; therefore, since the ADC has
4096 bins, the threshold in ADC bins must be
Thresh. =

Raw thresh.
× 0.175 × 4096,
128

(4.18)

where the raw threshold is the value determined from the two settings described above. The raw threshold is divided by 128 because that is the maximum value of a 7 bit unsigned integer.
Any time the amplified, shaped input signal rises above the threshold value, a trigger signal is generated by the AGET for the corresponding channel. The trigger signal is a square wave with a height of
48 ADC bins and a width of approximately2 235 ns. The trigger signal is either on or off: no trigger pileup
can occur at the channel level. If the signal is still above threshold at the end of the trigger pulse, another
trigger pulse will fire. This repeats until the signal falls below threshold.
2 This value is not well-defined, but is instead within a small tolerance due to the way the chip is designed [3].

66

4.1.3.2 Generating CoBo-level multiplicity signals
Before discussing the multiplicity signals generated by the CoBo, a brief detour is required to understand
the clocks used by the CoBo. The CoBo has a master clock frequency of 100 MHz from which several
other clock domains are derived. The two domains relevant to the analysis are the write clock and the
read clock:
• The write clock controls the sampling frequency in the detector, and its frequency is a configurable
fraction of the master clock. In the 46 Ar experiment, the write clock frequency was 12.5 MHz. When
the description of the analysis refers to a generic “clock frequency,” this is the frequency meant.
• The read clock controls the readout of the sampled values from the SCA into the ADC, and consequently the readout of digitized samples from the ADC to the CoBo. It has a fixed frequency of
25 MHz.
Between triggered events, each AGET sums the trigger signals from its 64 physics channels and continuously outputs this trigger multiplicity signal to the ADC. The CoBo sums the digitized multiplicity
signals from all of the AGETs it controls, and it integrates this signal over a sliding time window. The
multiplicity window has a configurable width which is defined in units of read clock ticks (or, units of
40 ns). This width was set to 300 ticks, or 12 µs, during the 46 Ar experiment. The integrated signal is then
compared to a threshold value, which was set to 20 000 during most of the experiment. If the integrated
multiplicity signal rises above this threshold, then the CoBo generates a trigger request.
An important caveat when simulating this process is that the CoBo samples the trigger signals from
each AGET using the read clock frequency, or every 40 ns. For convenience, the simulated trigger signals were sampled at the same 12.5 MHz frequency as the simulated events. Therefore, the simulated
integrated multiplicities must be multiplied a correction factor of
25 MHz
νr
=
=2
νw 12.5 MHz

(4.19)

since the CoBo would sample the signal twice as many times as the simulation during the multiplicity
window.

67

4.1.3.3 Generating global triggers
In the full GET electronics system, the channel-level trigger signals are collected by the MuTAnT module
and combined using a configurable algorithm to generate a global trigger for the detector. However, the
MuTAnT module was not yet available at the time of the

46

Ar(p, p ) measurement, so the AT-TPC was

triggered with an external system instead.
To simulate this external trigger, the CoBo-level trigger signals are each compared to their individual
integrated multiplicity thresholds. This produces a set of 10 Boolean values which are then combined
with a logical OR to generate the global trigger. Therefore, a global trigger is generated if any individual
CoBo rises above its multiplicity threshold during an event.

4.1.4 Cleaning and minimization
After testing whether the simulated events triggered the electronics, the events are processed by the
cleaning and fitting functions described above in Sections 3.3 and 3.5. The results of the fits are stored
along with the values of the parameters that were used to simulate the track.

4.1.5 Post-processing and cuts
After generating and fitting the tracks, the first step of the analysis process is to apply a cut on the values
of the objective function χ2 . Cuts are applied independently on the two components χ2pos and χ2en (see
Section 3.4.3 for definitions). The χ2 cuts also exclude events that cause the Monte Carlo fitter to return
an error. This happens, for example, if there are too few data points in the event.
A second, independent cut is applied on the results of the trigger step of the simulation. This is possible because, unlike for experimental data, the simulation code can still perform the Monte Carlo fit on
an event even if the event fails to fire the simulated trigger. This enables us to independently assess the
acceptances of the detector and the fitting software.
These two cuts divide the data set into six categories based on whether or not the detector was triggered and whether the fit was good, the fit was bad, or the fit failed with an error.

68

Pad threshold
Set

Runs

Pads excluded

Mult. thresh.

MSB

LSB

Notes

0
1
2
3
4
5

53–76
77
78–102
103–186
188–217
218–end

436
436
4247
436
665
665

15 000
15 000
20 000
20 000
20 000
20 000

1
1
1
1
1
1

6
2
2
2
2
2

CoBos 4 & 9 excluded

Table 4.3: Trigger conditions used in the experiment. As noted in the text, the unusually large number of
excluded pads in set 2 was due to a user error.

4.2 Simulated acceptance
Using the methods described above, a set of simulations was performed to estimate the acceptance of
the AT-TPC’s trigger setup and the efficiency of the fitting process. Because the trigger conditions were
adjusted several times during the experiment, a separate simulation was performed for each set of trigger
conditions so that the results could later be used to correct the recorded data for the efficiency. The
trigger conditions are listed in Table 4.3. The conditions for run set 1 were skipped since that set only
includes one run. For each of the remaining sets of runs, 106 events were simulated using a uniform
parameter distribution as described in Section 4.1.1.1.
One of the most significant differences between the run sets was the size and shape of the trigger exclusion and low-gain regions on the pad plane. These regions are shown in Fig. 4.4. The trigger exclusion
regions were defined by accumulating a large number of beam-only events and selecting the pads that
recorded the largest signals. The large size of the region for set 2 was due to a user error, and the wedge
shape of the region for set 5 was due to the exclusion of all of CoBos 4 and 9 from the trigger in an effort
to reduce the trigger rate.
The distributions of the accepted and rejected events are shown as functions of the track parameters
in Figs. 4.7 to 4.11. The trigger conditions appear to reject most events with a proton energy greater than
approximately 4 MeV or a center-of-mass scattering angle greater than approximately 80°. Interestingly,
the acceptance is not isotropic with respect to the azimuthal angle, but instead shows some fluctuations.
This happens because the azimuthal symmetry of the detector is broken when the detector is tilted and
placed inside the magnetic field.
In contrast with the χ2 distributions shown in Figs. 4.5 and 4.6, the acceptances shown in Figs. 4.7

69

Figure 4.4: Trigger exclusion and low-gain regions used in the simulation and in the experiment. Only
pads in the light and dark gray regions contribute to the trigger.

70

Count

4000
3000
2000
1000
0

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
pos

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
en

Count

(a) Set 0

2000
1500
1000
500
0

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
pos

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
en

Count

(b) Set 2

4000
3000
2000
1000
0

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
pos

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
en

(c) Set 3

Figure 4.5: Distributions of the χ2 values for the simulated events in the first three sets of runs. Events
that fell outside of the shaded region between the two vertical black lines on the histogram were thrown
away. The events that were kept had χ2pos < 1.0 and χ2en < 1.0. This plot does not include events for which
the fitter returned an error.

71

Count

4000
3000
2000
1000
0

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
pos

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
en

(a) Set 4

Count

3000
2000
1000
0

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
pos

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
en

(b) Set 5

Figure 4.6: Distributions of χ2 values for the final two sets of runs. See the caption of Fig. 4.5 for details.

72

6000
4000
2000
0

20

0

x [mm]

20

40

25000
20000
15000
10000
5000
0

0.0 2.5 5.0 7.5 10.0 12.5 15.0

Proton energy [MeV]

6000

Count

Count

Count
40

40

20

0

y [mm]

20

0

200

400

600

z [mm]

4000
2000
0

40

Count

20000
15000
10000
5000
0

Count

Count

20000
15000
10000
5000
0

800 1000

No trigger, fit failed
No trigger, bad fit
No trigger, good fit

10000
8000
6000
4000
2000
0

0

60 120 180 240 300 360

Azimuthal angle [deg]

30

60

90

120 150 180

CM angle [deg]

Trigger, fit failed
Trigger, bad fit
Trigger, good fit

Figure 4.7: Simulated acceptance of the trigger and analysis as a function of the track parameters for trigger parameter set 0. The trigger and analysis cuts were applied independently to the simulated events,
leading to six classes of events. Events are grouped into “trigger” and “no trigger” groups depending on
whether they would have triggered the electronics in a real experiment. These groups are then subdivided into “good fit” and “bad fit” groups depending on whether they pass the χ2 cuts. Events in the “fit
failed” groups caused an error in the Monte Carlo algorithm. Only events in the “trigger, good fit” group
would have been recorded and correctly reconstructed in the experiment. The six groups are stacked in
these histograms, so the contour of the top of each plot reflects the shape of the total distribution.

73

6000
4000
2000
0

20

0

x [mm]

20

40

0.0 2.5 5.0 7.5 10.0 12.5 15.0

Proton energy [MeV]

6000

Count

Count

Count
40

25000
20000
15000
10000
5000
0

40

20

0

y [mm]

20

0

200

400

600

z [mm]

4000
2000
0

40

Count

20000
15000
10000
5000
0

Count

Count

20000
15000
10000
5000
0

800 1000

No trigger, fit failed
No trigger, bad fit
No trigger, good fit

8000
6000
4000
2000
0

0

60 120 180 240 300 360

Azimuthal angle [deg]

30

60

90

120 150 180

CM angle [deg]

Trigger, fit failed
Trigger, bad fit
Trigger, good fit

Figure 4.8: Simulated acceptance of the trigger and analysis for trigger parameter set 2. The large trigger
exclusion region used in these runs prevented events outside of a very small angular and energy domain
from being accepted.

74

6000
4000
2000
0

Count

Count
20

0

x [mm]

20

40

25000
20000
15000
10000
5000
0

40

20

0

y [mm]

20

0

200

400

600

z [mm]

0.0 2.5 5.0 7.5 10.0 12.5 15.0

Proton energy [MeV]

Count

40

40

Count

20000
15000
10000
5000
0

Count

Count

20000
15000
10000
5000
0

800 1000

No trigger, fit failed
No trigger, bad fit
No trigger, good fit

6000
4000
2000
0

10000
8000
6000
4000
2000
0

0

60 120 180 240 300 360

Azimuthal angle [deg]

30

60

90

120 150 180

CM angle [deg]

Trigger, fit failed
Trigger, bad fit
Trigger, good fit

Figure 4.9: Simulated acceptance of the trigger and analysis for trigger parameter set 3.

75

6000
4000
2000
0

Count

Count
20

0

x [mm]

20

40

40

20

0

y [mm]

20

0

200

400

600

z [mm]

0.0 2.5 5.0 7.5 10.0 12.5 15.0

Proton energy [MeV]

Count

40

20000
15000
10000
5000
0

40

Count

20000
15000
10000
5000
0

Count

Count

20000
15000
10000
5000
0

800 1000

No trigger, fit failed
No trigger, bad fit
No trigger, good fit

6000
4000
2000
0

8000
6000
4000
2000
0

0

60 120 180 240 300 360

Azimuthal angle [deg]

30

60

90

120 150 180

CM angle [deg]

Trigger, fit failed
Trigger, bad fit
Trigger, good fit

Figure 4.10: Simulated acceptance of the trigger and analysis for trigger parameter set 4.

76

6000
4000
2000
0

20

0

x [mm]

20

40

0.0 2.5 5.0 7.5 10.0 12.5 15.0

Proton energy [MeV]

6000

Count

Count

Count
40

25000
20000
15000
10000
5000
0

40

20

0

y [mm]

20

0

200

400

600

z [mm]

4000
2000
0

40

Count

20000
15000
10000
5000
0

Count

Count

20000
15000
10000
5000
0

800 1000

No trigger, fit failed
No trigger, bad fit
No trigger, good fit

8000
6000
4000
2000
0

0

60 120 180 240 300 360

Azimuthal angle [deg]

30

60

90

120 150 180

CM angle [deg]

Trigger, fit failed
Trigger, bad fit
Trigger, good fit

Figure 4.11: Simulated acceptance of the trigger and analysis for trigger parameter set 5. The acceptance
was reduced for these runs due to the exclusion of CoBos 4 and 9 from the trigger. This is most visible in
the distribution with respect to the azimuthal angle.

77

to 4.11 also include events that were not successfully processed by the Monte Carlo code. The majority of
these events, which are shown in the “trigger, fit failed” and “no trigger, fit failed” groups, failed because
they had a very small number of data points remaining after cleaning. Specifically, the fitting algorithm
rejects any event with fewer than 50 data points without attempting to fit it as these events tend to be
short tracks that are fit poorly. These rejected events are therefore not shown in the χ2 distributions since
no χ2 values were calculated for them.
This simulated data set was used to calculate the efficiency of the detector as a function of centerof-mass scattering angle and

46

Ar vertex energy. To this end, the data were first binned in those two

dimensions. The bins were assigned based on the known correct values of the scattering angle and vertex
energy used to simulate the tracks—rather than the values reconstructed by the Monte Carlo fit—in order
to remove any uncertainties arising from the fit itself. The efficiency in a given 2D bin i was defined as

i

=

Ai
,
Ni

(4.20)

where A i is the number of accepted events in the bin and Ni is the total number of events simulated in
the bin. An “accepted” event can be defined as one that triggers the detector, passes the χ2 cut, or both.
Two-dimensional histograms of the efficiencies are shown in Figs. 4.12 and 4.13.
Statistically, the acceptance or rejection of an event within a bin can be considered to be a Bernoulli
random process: the event is accepted with a fixed probability , which is equal to the efficiency, or rejected with probability 1 − . This implies that the number of accepted counts A i in each bin should be
binomially distributed with variance Ni (1 − ), where Ni is, again, the number of events simulated in
that bin. Since the true value of

is unknown, we can approximate the variance with the calculated

i.

This leads to the following formula for the uncertainty of the efficiency:
σ

,i

=

σA
1
=
Ni
Ni

Ai 1 −

Ai
.
Ni

(4.21)

Figures 4.14 to 4.18 show the calculated efficiencies with these error bars as a function of vertex energy
for a variety of scattering angles.
The plots shown in this section suggest that the limiting factor for the acceptance in this experiment
was the trigger. The low acceptance of the trigger dominated the acceptance of the track reconstruction
and fit, especially at forward angles. There are a few factors that likely contributed to this effect:

78

Analysis

Total

1.0

4

0.8

3

0.6

2

0.4

1

0.2

0

Efficiency

46Ar vertex energy [MeV/u]

Trigger

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

0.0

CM angle [deg]

Trigger

Analysis

Total

1.0

4

0.8

3

0.6

2

0.4

1

0.2

0

Efficiency

46Ar vertex energy [MeV/u]

(a) Set 0

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

0.0

CM angle [deg]

Trigger

Analysis

Total

1.0

4

0.8

3

0.6

2

0.4

1

0.2

0

Efficiency

46Ar vertex energy [MeV/u]

(b) Set 2

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

0.0

CM angle [deg]

(c) Set 3

Figure 4.12: Efficiencies of the AT-TPC trigger, the analysis χ2 cuts, and the total of the two for the first
three sets of parameters, plotted as functions of the 46 Ar energy at the reaction vertex and the scattering
angle in the center of mass frame. The bins used were 5° wide in the angular dimension and 50 keV wide
in the energy dimension.

79

Analysis

Total

1.0

4

0.8

3

0.6

2

0.4

1

0.2

0

Efficiency

46Ar vertex energy [MeV/u]

Trigger

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

0.0

CM angle [deg]

Trigger

Analysis

Total

1.0

4

0.8

3

0.6

2

0.4

1

0.2

0

Efficiency

46Ar vertex energy [MeV/u]

(a) Set 4

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

CM angle [deg]

40

60

80 100 120

0.0

CM angle [deg]

(b) Set 5

Figure 4.13: Efficiencies of the AT-TPC trigger, the analysis χ2 cuts, and the total of the two for the last
two sets of parameters. See the caption of Fig. 4.12 for details.

80

Efficiency

1.00

CM = 40°

CM = 45°

CM = 50°

CM = 55°

CM = 60°

CM = 65°

CM = 70°

CM = 75°

CM = 80°

CM = 85°

CM = 90°

CM = 95°

0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00
0

1

2

3

4

46Ar vertex energy [MeV/u]

0

1

2

3

4

46Ar vertex energy [MeV/u]

Trigger

Analysis

0

1

2

3

4

46Ar vertex energy [MeV/u]

Total

Figure 4.14: Efficiencies for parameter set 0 as a function of 46 Ar vertex energy at several scattering angles.
The error bars shown are ±σ ,i as defined in Eq. (4.21). The bins used in this plot are the same as those
used in the 2D histograms in Figs. 4.12 and 4.13.

81

Efficiency
Efficiency
Efficiency
Efficiency

0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
0.0
0

1

CM = 40°

CM = 45°

CM = 50°

CM = 55°

CM = 60°

CM = 65°

CM = 70°

CM = 75°

CM = 80°

CM = 85°

CM = 90°

CM = 95°

2

3

4

46Ar vertex energy [MeV/u]

0

1

2

3

4

46Ar vertex energy [MeV/u]

Trigger

Analysis

0

1

2

3

4

46Ar vertex energy [MeV/u]

Total

Figure 4.15: Efficiencies for parameter set 2 as a function of 46 Ar vertex energy at several scattering angles.

82

Efficiency

1.00

CM = 40°

CM = 45°

CM = 50°

CM = 55°

CM = 60°

CM = 65°

CM = 70°

CM = 75°

CM = 80°

CM = 85°

CM = 90°

CM = 95°

0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00
0

1

2

3

4

46Ar vertex energy [MeV/u]

0

1

2

3

4

46Ar vertex energy [MeV/u]

Trigger

Analysis

0

1

2

3

4

46Ar vertex energy [MeV/u]

Total

Figure 4.16: Efficiencies for parameter set 3 as a function of 46 Ar vertex energy at several scattering angles.

83

Efficiency

1.00

CM = 40°

CM = 45°

CM = 50°

CM = 55°

CM = 60°

CM = 65°

CM = 70°

CM = 75°

CM = 80°

CM = 85°

CM = 90°

CM = 95°

0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00

Efficiency

1.00
0.75
0.50
0.25
0.00
0

1

2

3

4

46Ar vertex energy [MeV/u]

0

1

2

3

4

46Ar vertex energy [MeV/u]

Trigger

Analysis

0

1

2

3

4

46Ar vertex energy [MeV/u]

Total

Figure 4.17: Efficiencies for parameter set 4 as a function of 46 Ar vertex energy at several scattering angles.

84

Efficiency

0.8

CM = 40°

CM = 45°

CM = 50°

CM = 55°

CM = 60°

CM = 65°

CM = 70°

CM = 75°

CM = 80°

CM = 85°

CM = 90°

CM = 95°

0.6
0.4
0.2
0.0

Efficiency

0.8
0.6
0.4
0.2
0.0

Efficiency

0.8
0.6
0.4
0.2
0.0

Efficiency

0.8
0.6
0.4
0.2
0.0
0

1

2

3

4

46Ar vertex energy [MeV/u]

0

1

2

3

4

46Ar vertex energy [MeV/u]

Trigger

Analysis

0

1

2

3

4

46Ar vertex energy [MeV/u]

Total

Figure 4.18: Efficiencies for parameter set 5 as a function of 46 Ar vertex energy at several scattering angles.

85

1. Higher-energy protons do not generally stop inside the active volume of the detector. Instead, they
spiral through the volume at a nearly constant energy, far from the peak of the Bragg curve, and
free fewer electrons per unit track length. This creates smaller signals in the electronics that are
less likely to rise above threshold, preventing individual channels from triggering.
2. Events at very forward angles illuminate fewer pads on the sensor plane, leading to lower multiplicities. In principle, tilting the detector should address this since it spreads out the projection of
the beam—and therefore the projection of forward events—over more pads; however, the particular combination of drift velocity, tilt angle, and magnetic field in this experiment combined to
produce a Lorentz angle that almost perfectly canceled the tilt.
3. Since the MuTAnT module was not available during the experiment, the trigger used was not a
true global multiplicity trigger, but instead a set of CoBo-level multiplicity triggers combined using a logical OR operation. Thus, events that activated few pads in many CoBos were cut if the
multiplicity in each CoBo was below threshold, even if the event had a high global multiplicity.
Many of these effects can be addressed in future experiments by using the MuTAnT as a trigger and by
carefully choosing the combination of gas pressure and fields.

4.3 Simulated resolution
The simulated data set described in Section 4.2 can also be used to establish a benchmark for the expected resolution of the detector for 46 Ar(p, p ) scattering events. This was done by finding the difference
between the reconstructed track parameters and the known correct parameters that were used to simulate each event. The distributions of this difference for the six track parameters are shown in Fig. 4.19.
Table 4.4 shows the output of fitting functions to these distributions of deviations. The distributions
for ∆x 0 , ∆y 0 , ∆E pos , and E kin − E pos were fit using a Laplace distribution, which was defined as
f (x, µ, b, A) =

|x − µ|
A
exp −
.
2b
b

(4.22)

In the nomenclature used in Table 4.4, µ serves as the centroid for this distribution, b is the width parameter, and A is the normalization parameter. The remaining distributions were fit with a Gaussian

86

Count

Count

1e3
6
4
2
0

2.0
1.5
1.0
0.5
0.0

40

20

0

x0 [mm]

20

0.8
0.6
0.4
0.2
0.0

40

1e3

100 50

1e4

6
4
2

40

20

0

y0 [mm]

20

40

0

1e3

0

50 100

E0 [keV]

5
4
3
2
1
0

1e3

40

40

20

0

20

40

0.5 0.0

0.5

1.0

z0 [mm]

1e3

20

0

0

[deg]

20

40

3
2
1
0

1.0

0

[deg]

Figure 4.19: Simulated resolution of the AT-TPC. These six plots show the error in the reconstruction of
each of the six track parameters, defined as the reconstructed value minus the correct value. The results
of the fits to these distributions are given in Table 4.4.

Variable

Distribution

Centroid

Width

FWHM

Units

Normalization

∆x 0
∆y 0
∆z 0
∆E 0
∆φ0
∆θ0
∆E kin
∆E pos
E kin − E pos

Laplace
Laplace
Gaussian
Gaussian
Gaussian
Gaussian
Gaussian
Laplace
Laplace

−0.210
−0.245
6.19
1.18
2.06 × 10−2
3.05 × 10−2
−7.25
20.4
−48.0

5.07
4.30
4.81
14.3
5.00
0.144
244
18.0
229

7.03
5.96
11.3
33.7
11.8
0.339
575
25.0
317

mm
mm
mm
keV
deg
deg
keV/u
keV/u
keV/u

6.72 × 104
6.32 × 104
6.32 × 104
6.21 × 104
6.39 × 104
1.20 × 103
6.64 × 105
1.36 × 105
7.14 × 105

Table 4.4: Results of fitting the distributions of parameter and energy deviations from Figs. 4.19 and 4.21.
The Laplace and Gaussian distributions are defined in Eqs. (4.22) and (4.23), respectively. The units listed
apply to both the centroid and width parameters. The full width at half-maximum (FWHM) was defined
as 2σ 2 log(2) for the Gaussian distributions and 2b log(2) for the Laplace distributions.

87

y

∆y 0

x

∆x 0

Figure 4.20: Illustration of relationship between uncertainty in vertex location and error in azimuthal angle. The reconstructed vertex is marked with an orange square, while the actual vertex location is marked
with a blue circle. The red curve represents the proton track. The tangent lines at these two locations
have very different slopes.
distribution
f (x, µ, σ, A) =

A
2σ2 π

exp −

(x − µ)2
,
2σ2

(4.23)

where µ is the centroid, σ is the width parameter, and A is the normalization parameter. These results
show a FWHM resolution of 33.7 keV for the proton energy and 0.339° for the scattering angle in the
laboratory frame. The z position of the proton was resolved to within 11.3 mm at FWHM.
One unusual aspect of these results is the poor resolution found for the azimuthal angle φ0 . The
error in φ0 is caused by the uncertainty in the extrapolation from the first data point back to the reaction
vertex. The uncertainty in the vertex location leads to a correlated uncertainty in the azimuthal angle.
This is shown by the difference in the tangent lines at the actual vertex and the reconstructed vertex in
Fig. 4.20, which corresponds to ∆φ0 .
The parameters of the simulated tracks can also be used to estimate how precisely the 46 Ar vertex energy can be reconstructed using the same techniques as with data (see Section 3.7). The kinematic and
positional vertex energies can be calculated for both the reconstructed parameters and the parameters
used to simulate the tracks, providing an estimate of the resolution. The results of this comparison are
plotted in Fig. 4.21, which shows distributions of the deviation in kinematic vertex energy and positional
vertex energy as well as the difference between the two. The kinematic vertex energy was reconstructed

88

1e3

1e3

1.0

Count

1e3
1.5

3
1.0

2

0.5

0.5

1
0.0

0
1000 500

0

500 1000

Ekin [keV/u]

0.0
150 100 50 0 50 100 150

1000 500

Epos [keV/u]

0

500 1000

Ekin Epos [keV/u]

Count

Figure 4.21: Simulated resolution of the two energy reconstruction techniques from Section 3.7. The
widths of the distributions are given in Table 4.4. ∆E kin and ∆E pos are defined as the difference between
the energies as calculated from the reconstructed parameters and the known correct parameters. For
the rightmost plot of E kin − E pos , both values were found from the reconstructed parameters. The difference between the two energies for the known correct parameters is not shown since those two values are
identical. Note that the horizontal scale is different on the center plot.

10000
8000
6000
4000
2000
0

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
pos

0.0 0.5 1.0 1.5 2.0 2.5 3.0
2
en

Figure 4.22: Distributions of the objective function components for simulated events drawn from the R
matrix calculation results. Each distribution was cut at 1.0, discarding all events with a larger value of χ2 .
This is the same cut that was used for the flat distribution (see Figs. 4.5 and 4.6).
to a resolution of 575 keV/u, the positional vertex energy was resolved to within 42.4 keV/u, and the difference between the two had a resolution of 539 keV/u, all at FWHM.

4.4 R matrix simulation
A second simulation was performed to assess the AT-TPC’s ability to resolve the resonances that were expected in the 46 Ar(p, p ) measurement. This was done by simulating approximately 5.5 × 105 events following a distribution derived from an R matrix calculation as described in Section 4.1.1. Of these events,
approximately 3.4 × 105 passed the simulated trigger cuts and the χ2 cuts shown in Fig. 4.22.

89

5000

Count (normalized w.r.t. efficiency)

46Ar vertex energy [MeV/u]

4.200
3.875

4000
3000

3.000
2.696

2000
1000

2.000

40

50

60

70

80

90

CM angle [deg]

100

110

120

0

Figure 4.23: Two-dimensional histogram of results from the R matrix–based simulation. The two expected resonance locations are marked on the vertical axis at 2.696 MeV/u and 3.875 MeV/u. The two
individual bright bins at 80° and 85° are artifacts from dividing by the very small number of counts in
those bins.
After applying the χ2 cut, the remaining events were binned based on their reconstructed positional
46

Ar vertex energies and center-of-mass scattering angles. The reconstructed values were used instead of

the known correct values since the purpose of this particular simulation was to see how well the detector
and fit could measure the resonances. The binned data were then normalized with respect to the total
efficiency calculated earlier in Section 4.2.
Figure 4.23 shows a two-dimensional histogram of the results of the simulation. The resonances corresponding to the ground state and first excited state are visible as faint horizontal depressions in the
color scale at energies of 2.696 MeV/u and 3.875 MeV/u, respectively. These energies differ by a small
amount from the resonance energies listed in Table 4.1 since those correspond to the incoming proton
energies in the center-of-mass frame while the simulated energies in Fig. 4.23 correspond to the 46 Ar energies in the laboratory frame. The two values are related by Eqs. (4.12) and (4.16) from Pages 63 and 64.
Figure 4.24 shows excitation functions from this simulation taken at a number of different centerof-mass angles. The resonances are particularly visible in the excitation functions for center-of-mass
angles of 55° and 60°. This is partially because these are some of the largest angles where the efficiency

90

of the trigger is still large enough to accept many events. Larger angles correspond to smaller impact
parameters and more-central collisions, and therefore should show larger resonances as seen in the R
matrix calculation results in Fig. 4.2. However, when the efficiency becomes too small, the low statistics
prevent a conclusive determination of whether or not the resonance is observed. The statistics at higher
center-of-mass angles are naturally decreasing since the Rutherford cross section approaches zero as
θCM approaches 180°.

91

Count (normalized)

6000

CM = 40°

4000

2000

2000

2000

1000

0

0

0

2000

2000

Count (normalized)

CM = 60°

1500

500

500

0

0

1000
750
500
250
0

0
CM = 75°

6000

1000

4000

500

2000

0

0
CM = 90°

1000

2000

500

1000

CM = 95°

600
400
200

0
2.0 2.5 3.0 3.5 4.0

CM = 80°

800

3000

46Ar vertex energy [MeV/u]

8000

1500

CM = 85°

1500

CM = 65°

1000

1000

1000

CM = 70°

Count (normalized)

CM = 50°

3000

4000

CM = 55°

Count (normalized)

CM = 45°

2.0 2.5 3.0 3.5 4.0

46Ar vertex energy [MeV/u]

2.0 2.5 3.0 3.5 4.0

46Ar vertex energy [MeV/u]

Figure 4.24: Excitation functions derived from the R matrix simulation. The expected locations of the two
resonances are shown with vertical orange lines.

92

CHAPTER 5
EXPERIMENT

This chapter describes an experiment that was performed at the National Superconducting Cyclotron
Laboratory (NSCL) to measure the reaction 46 Ar(p, p ). The experiment, which was assigned the NSCL experiment number e15503b ran from the morning of September 9, 2015, through the evening of September 20, 2015. The description begins with an overview of the beam production methods and the detector
setup used to measure the reaction. Following that is a description of the parts of the analysis process
that are particular to this experiment; the general analysis process was described in Chapter 3. Finally,
the experimental results are shown and discussed.

5.1 Beam production
The beam of

46

Ar used in this experiment was produced by the Coupled Cyclotron Facility [34] at the

National Superconducting Cyclotron Laboratory (NSCL) at Michigan State University. A diagram of this
facility’s beam line is shown in Fig. 5.1. Beam production began with a primary beam of 48 Ca produced by
a superconducting ion source and accelerated to a final energy of 140 MeV/u by the coupled cyclotrons.
The primary beam was then impinged upon a beryllium target to produce a mixture of nuclei by the
process of projectile fragmentation. These nuclei were separated in the A1900 fragment separator (see
Fig. 5.2), a device which uses four large dipole magnets to separate particles based on their magnetic
rigidity [36]. Particles with a different rigidity from the desired

46

Ar fragments follow different paths

through the separator and are filtered out of the beam using adjustable slits.
The 46 Ar beam exited the A1900 at a much higher energy than what was needed to see the resonances
in this experiment. However, the primary beam cannot be accelerated to a lower initial energy since at
very low energies, fragmentation is not possible. This issue was addressed by stopping the high-energy
beam particles and then re-accelerating them to a lower final energy. To accomplish this, the beam was
sent into a helium-filled [53] linear gas cell where it was stopped. Electric fields were then used to extract
the stopped ions from the gas cell and transport them to an electron beam ion source (EBIS), where their
charge states were increased from 1+ to 17+ to make acceleration possible. Finally, the ions were injected
into ReA3, a linear accelerator (linac), and re-accelerated to a final energy of 4.6 MeV/u.

93

Figure 5.1: A diagram of the NSCL beamline, with some regions of the facility labeled. The beam starts at
the ion sources (label 1) and ends at the AT-TPC (label 14). (Reproduced from http://nscl.msu.edu/
public/virtual-tour.html.)

91

D.J. Morrissey et al. / Nucl. Instr. and Meth. in Phys. Res. B 204 (2003) 90–96
ion sources

K1200

RT-ECR

20 ft

SC-ECR

A1900

image 1
image 2

stripping
foil

production
target

10 m

K500

coupling
line

focal plane
image 3

Fig. 1. A schematic diagram of the coupled cyclotron facility and the operation of the A1900 fragment separator, see the text.

Figure 5.2: The A1900 fragment separator and the coupled cyclotrons. (Reproduced from [36].)

Fig. 1 shows a schematic layout of the A1900
2000 and the commissioning studies of the accelseparator connected to the K500 and K1200 cyerators and fragment separator started in 2001.
clotrons. Intense beams of stable ions from the
The general features of the new equipment are
electron cyclotron-resonance ion sources are indescribed below and the results of initial studies of
jected into the center of the K500 cyclotron and
the parameters of the separator are also presented.
94
accelerated to an energy sufficient to strip most of
Several experiments with exotic secondary beams
the remaining electrons in a foil positioned near
have already been completed.
the center of the K1200 cyclotron. The highly
The heart of any nuclear physics program with

40

Event ID

30
20
10
0
0

1

2

3

Time [s]

4

5

6

Figure 5.3: Timing of events recorded by the GET electronics. Particles arrived in bunches at a rate of
approximately 2 Hz. Note that this only shows events that triggered the electronics; in reality, many more
nuclei arrived in each bunch.
This re-acceleration process is very useful as it allows rare isotope beams to be produced at much
lower energies than was previously possible; however, the re-accelerated beams do have a few unique
characteristics that must be accounted for. The primary complication introduced is the time structure of
the re-accelerated beam. Particles are extracted from the EBIS with a period on the order of milliseconds,
creating large, widely spaced bunches of particles (see Fig. 5.3). Additionally, the first stage of the ReA3
linac was operated in a pulsed mode to address electrical power considerations, further contributing to
bunching. This resulted in a very high instantaneous beam rate with a low duty cycle, causing a significant degree of pileup in the AT-TPC, as shown in Fig. 5.4. These high-multiplicity events are very difficult
to analyze, but fortunately they can be easily filtered out using data from an ion chamber upstream of
the AT-TPC, as described in Section 5.3.3 below.
In addition to bunching the beam, the duration of the charge breeding process allows some of the
beam nuclei to decay. While the 8.4(6) s half-life [38] of 46 Ar is very long compared to the time it takes
the particles to travel from the production target to the EBIS, the nuclei then spend on the order of milliseconds in the EBIS. During that time, a small but significant number of 46 Ar nuclei β-decay into 46 K,
a few of which further β-decay into 46 Ca with a half-life of 105(10) s [38]. These decay products have the
same mass as 46 Ar and some of them leave the EBIS with the same charge state, so they are transmitted
through the pre-separator and accelerated to the same final velocity by the linac. Fortunately, it is possible to separate the 46 Ar from its decay products using the signals from the ion chamber, which will also

95

IC signal amplitude

3000
2500
2000
1500
1000
500
0
0

100

200

300

ADC time bucket

400

500

Figure 5.4: A signal from the ion chamber (see Section 5.2.2) from an event where four particles entered
the detector.
be described in Section 5.3.3 below.

5.2 Detection setup
The beam of

46

Ar was measured using the AT-TPC, which was described in detail in Chapter 2. This

section will cover the way the detector was set up in this experiment.

5.2.1 AT-TPC setup
The AT-TPC was mounted inside its solenoid in the ReA3 experimental area at the NSCL and connected
to the ReA3 beam line. The jack on the detector’s support carriage was used to raise the downstream end
of the detector until its central axis was tilted upward by an angle of 6.2° with respect to the central axis
of the solenoid. The GET electronics were then connected to the TPC.
The active volume of the detector was filled with isobutane gas at a pressure of 19.2 torr and the
shielding volume was filled with nitrogen gas at the same pressure. Both gases were circulated through
the detector and continuously replaced at a rate of 83 cm3 /min to maintain a low level of contaminants
in the gas.
A potential of −9.0 kV was applied to the AT-TPC’s cathode with respect to its grounded anode, creating an electric field of approximately 9.5 kV/m in the active volume. The solenoid was energized and
produced a magnetic field of 0.997 T as measured by a Hall probe mounted on the axis of the magnet and

96

flush with the upstream end of the bore. According to a field map of the magnet, this corresponds to a
magnetic field of 1.68 T in the center of the bore where the TPC is located. The solenoid’s power supply
was disconnected after the magnet was energized, leaving the magnet running in persistent mode.
To amplify the ionization electrons, a potential of −450. V was applied to the Micromegas mesh with
respect to the grounded pad plane. In later runs, this was reduced to −446. V to reduce sparking between
the mesh and the pads. This high potential difference was required to detect the relatively small amount
of charge produced by the scattered protons even though it caused a relatively high current of around
8 µA to flow between the pads and the mesh when the heavily ionizing 46 Ar nuclei passed through the
detector. A potential of −500. V was applied to the field cage ring nearest to the Micromegas to ensure
the linearity of the electric field between the end of the field cage and the negatively biased Micromegas
mesh.

5.2.2 Ion chamber and secondary DAQ
A small, cylindrical ion chamber was constructed and inserted into the beam line just upstream of the
AT-TPC. The chamber is 5 cm long and contains 5 thin foils mounted at equal intervals along the axial direction. The foils have a diameter of 1 in and consist of approximately 400 µg/cm2 of aluminum
evaporated on a plastic substrate. The outermost two foils are grounded and serve as entrance and exit
windows, while the inner three foils establish equipotentials for the electric field. The electron drift direction is parallel to the beam, and signals are read out from the center foil, which is kept at a potential
of 100 V relative to the outermost foils. Like the active volume of the AT-TPC, the ion chamber was filled
with isobutane gas at a pressure of 19.2 torr.
This ion chamber served two purposes. First, the energy loss data from the ion chamber was used
to help separate the desired beam particles from contaminants. Second, the timing signal from the ion
chamber was used to prevent the AT-TPC from triggering on events that were not correlated with a beam
particle and to establish an absolute, global time reference for each event that corresponded to when
a beam particle entered the active volume of the AT-TPC. The data from the ion chamber was digitized
using a 12-bit flash ADC VME module and recorded using a separate DAQ from the main system. This
auxiliary DAQ system also digitized and recorded the trigger signal output to the CoBos and a signal
derived from the Micromegas mesh.

97

5.2.3 Trigger
The final version of the MuTAnT board was not yet available when this experiment was run, so the GET
electronics had to be triggered externally. Each CoBo was configured to compare its total multiplicity
signal to a multiplicity threshold with an integration window of 300 cycles of the 25 MHz read clock, or
12 µs. Whenever the CoBo multiplicity exceeded the threshold, the CoBo would output a trigger signal
through a LEMO connector on its front panel. The CoBos also output a second signal through another
connector whenever they were busy processing an event. For a more thorough description of the CoBo
multiplicity trigger, see Section 4.1.3. The trigger thresholds are listed in Table 4.3 on Page 69.
To build a trigger, the multiplicity trigger signals from the 10 CoBos were collected and combined
with a logical OR operation in a VME module to produce a global multiplicity trigger signal. The busy
signals were combined in the same way to produce a global busy signal. These two global signals were
then combined to produce a signal that was true when any CoBo had a multiplicity trigger and no CoBo
was busy. More formally, for CoBo multiplicity signals M i and busy signals B i , the trigger logic derived
from the GET electronics was
TGET = (M 0 ∨ M 1 ∨ · · · ∨ M 9 ) ∧ (B 0 ∨ B 1 ∨ · · · ∨ B 9 ).

(5.1)

In addition to the signals from the GET electronics, a trigger signal was derived from the ion chamber
output using a discriminator module. This ion chamber trigger was then delayed and compared to a
coincidence window which was opened by the CoBo multiplicity trigger. If the delayed ion chamber
signal fell within this coincidence window, a master trigger signal was generated and distributed to the
10 CoBos.
The master trigger was based on the delayed ion chamber signal to ensure that the trigger timing
was always constant with respect to the arrival of the beam particle into the AT-TPC, rather than the
delayed time of arrival of the drifting electrons. The latter is highly dependent on the trajectory of the
charged particle emitted during the reaction, and therefore it does not provide a good timing reference.
This reference is essential to establish the relationship between the measured time of the signals and the
position along the drift direction of the electrons since this relationship is used to deduce the absolute
position of the tracks in this direction.
Finally, in addition to the master trigger above, the detector was also periodically triggered on a
downscaled version of the raw trigger signal from the ion chamber. This was done so that the detector
98

would record a few unreacted beam events every second for calibration purposes.

5.2.4 Excluding beam events
The GET electronics provide the ability to remove individual pads from the multiplicity trigger described
in Section 5.2.3 and to adjust the gain of the charge-sensitive preamplifier on each channel independently (see Fig. 4.3 on Page 66 for a schematic of the channel, where the preamplifier is labeled “CSA”).
These features are crucial as the pad-level trigger suppression allows unreacted beam events to be excluded from the trigger, and the channel-level preamplifier gain setting allows a low-gain region to be
established where the strongly ionizing beam particles are projected onto the sensor plane.
These regions were found by taking a short data run with no beam suppression. Approximately 100
events were then analyzed and the integrated charge deposited in each pad was accumulated for all
events into one array. An empirically chosen threshold was applied to select the set of pads that were
activated most often in this set of events. Since scattering events are rare compared to unreacted beam
events, this set of pads corresponds to the set of pads illuminated by the beam. These pads were then set
to a low gain and removed from the trigger. The exclusion regions used in the experiment were shown in
Fig. 4.4 on Page 70.

5.3 Analysis and results
The analysis process for AT-TPC data was described in general in Chapter 3. This section will cover the
specifics of the analysis done on the 46 Ar data, and then the results of the analysis will be shown.

5.3.1 Correlating VME and GET events
As described in Section 5.2.2, a separate VME DAQ was used to record signals from the ion chamber and
Micromegas mesh. The ion chamber signal is required to select 46 Ar events and exclude beam contaminants, so each event recorded by the VME DAQ must be paired with its corresponding event in the GET
DAQ. In principle, this can easily be done by comparing the event ID recorded in each system: the two
IDs should be perfectly correlated since the VME DAQ generated the trigger signals that were sent to the
GET DAQ. However, in practice, the two sets of event IDs tended to drift out of alignment over the course
of longer runs.
99

2.0

3072

1.5

GET signal sum

GET signal ampl.

4096

2048
1024
0

1e6

1.0
0.5
0.0

Mesh signal ampl.

IC signal ampl.

8192
4096
3072
2048
1024
0

6144
4096
2048
0

100 0

128

256

384

GET ADC time bucket

512

100 0

128

256

384

GET ADC time bucket

512

Figure 5.5: Plots of various signals in a spark-like event. The top left plot shows the individual signals in
the GET electronics, and the top right plot shows the sum of these signals. The bottom left plot contains
the signal from the ion chamber, and the bottom right plot shows the signal from the Micromegas mesh.
The gray regions indicate domains in which there is no data. The time buckets in the two DAQ systems
are of the same width, but are offset by 100 due to the trigger delay.

Uncalibrated y [mm]

150

CoBo 0
CoBo 1
CoBo 2
CoBo 3
CoBo 4
CoBo 5
CoBo 6
CoBo 7
CoBo 8
CoBo 9

100
50
0
50
100
150

100

0

100

Uncalibrated x [mm]

0

100

200

300

400

500

Uncalibrated z [GET ADC time buckets]

Figure 5.6: Plots of (x, y, z) positions corresponding to peak locations in a spark-like event. The majority
of the channels in the inner region of the pad plane fired at approximately the same time.

100

The alignment of the two systems was checked and corrected by looking at spark-like events. These
occurred periodically in every run, and were characterized by a large saturated pulse in the signal from
the mesh (the bottom right plot in Fig. 5.5) and large, simultaneous pulses in many channels of the GET
electronics (Fig. 5.6).
A small Python script was developed to analyze each spark-like event in the two DAQs and realign the
event IDs when one system became misaligned with the other. Spark-like events in the VME DAQ were
identified by selecting events whose mesh signal saturated the ADC (i.e., had a minimum value of 0) and
had a large, negative slope in a region more than one time bucket wide. The condition on the magnitude
of the slope removes pileup events from consideration, and the requirement that the region of large
slope be wider than one time bucket removes events where the signal suffered a discontinuity due to a
malfunction of the electronics. Any sets of sequential spark-like events were removed from consideration
since it would be impossible to determine whether the event in each DAQ was due to the same spark.
Once the spark-like events were identified in the VME DAQ, the GET event with the corresponding event ID was examined to determine if it contained a spark as well. Sparks in the GET events were
identified by selecting the region of time buckets where the mesh sparked and counting the number of
channels that attained an extreme value in that region. Extreme values were defined as values greater
than 4060 or less than 10 (for comparison, recall that the range of the 12-bit ADC is from 0 to 4095).
Channels that did not vary from their baseline by more than 10 time buckets were discarded. The event
was classified as a spark if more than 40 % of the remaining channels attained an extreme value. Outliers
in this statistic were identified by analyzing the two spark events following the current one if misalignment was detected: if the next two sparks were aligned, the current spark was labeled as an outlier and
ignored.
When the two systems were found to be out of alignment, the two GET events before and after the
event where the spark was expected were analyzed to try to find the correct alignment point. If this succeeded, an event ID offset was recorded for the remaining events in the run, and the events between the
last aligned spark and the current spark were marked as invalid since there is no way to determine where
in that range the misalignment began. If the alignment point could not be found, an error message was
printed and the run was deemed unrecoverable. The portion of each run between the last spark and the
end was marked as invalid because it would be impossible to tell if these events were misaligned.
Finally, note that this complex procedure will not be required in future experiments since the Mu-

101

TAnT board provides a global clock signal and generates the timestamps and event IDs used across the
entire system, thereby ensuring synchronization. As stated previously, though, the MuTAnT was unfortunately not yet available when the 46 Ar experiment was performed.

5.3.2 Merging GET data files
As mentioned in Section 2.5, the DAQ system features a distributed design that routes data from each
CoBo to an independent recording computer. While this improves throughput, it also means that the
experimental data for each run is spread across many separate files. The partial events from these files
had to be merged before they could be analyzed. The data format written by the GET system provides
two ways to merge events: by event ID and by timestamp.
The most reliable way to merge events is by timestamp since the timestamps should increment uniformly between CoBos even if not all CoBos are triggered on each event. However, merging by timestamp
was not possible for this experiment since the MuTAnT board was not available, and without the MuTAnT
providing a system-wide clock signal, the timestamps on each CoBo slowly diverged due to the natural
variance in the circuits’ internal clock frequencies.
Instead, the events were merged using the event ID. This method is generally reliable for two reasons.
First, the external trigger setup always distributed a trigger signal to all CoBos regardless of which CoBos
had a multiplicity trigger, so it should be impossible for any CoBo to not receive a trigger for a particular event. Second, the CoBos were configured to always output the fixed-pattern noise channels even if
there were no physics channels above threshold in a given event, so each data file will always contain at
least one data frame for every event even if that event did not contain any physics data for that CoBo.
That being said, during the analysis it became apparent that CoBo 0 occasionally failed to trigger, causing its event IDs to fall behind those produced by the rest of the CoBos. This was corrected by comparing
the intervals between successive events’ timestamps in each CoBo and inserting an offset where those
intervals were not approximately equal. This method works even though the CoBo clocks were not synchronized because the slight difference in clock periods between the boards is negligible on the small
scale of the interval between events.
The raw output was merged into one file per experimental run in a separate step before the rest of
the analysis was performed. The merging program read in partial events from the raw files produced by

102

8000

ADC bin

7000
6000
5000
4000
3000
2000

0

128

256

Time bucket

384

512

Figure 5.7: A raw ion chamber signal produced by an 46 Ar nucleus. The passage of the nucleus through
the chamber is indicated by the large negative peak around t = 160. The ADC used to record this signal was 12 bits, but it also had an overflow bit that made it effectively record 13 bits per sample for a
maximum digitized value of 8192.
the DAQ, grouped them by event number, and assembled full events from the fragments. The four fixedpattern noise (FPN) channels in each AGET were then averaged together and renormalized to produce a
set of AGET-level FPN signals with mean zero. These AGET-level FPN signals were then subtracted from
the physics channels controlled by each AGET to correct the data for the fixed-pattern noise. The FPN
data was then discarded. Finally, the data was compressed and written to an HDF5 file1 for later analysis.

5.3.3 Beam particle identification
Although the beam of 46 Ar delivered by ReA3 was about 90 % pure, there were a few contaminants that
needed to be removed from the data set. The main contaminants were identified as 46 K and 46 Ca (the
daughter and granddaughter nuclei of 46 Ar) as well as one or more isotopes of nickel from the beam line.
These were removed using the data from the ion chamber.
A raw signal from the ion chamber is shown in Fig. 5.7. The amplitude of the peak produced by the
beam particle is proportional to the amount of charge deposited by the particle as it passed through the
isobutane gas inside the chamber. This energy loss is, in turn, proportional to the square of the charge of
the beam particle according to Bethe’s formula. If we assume that all of the species present in the beam
have the same momentum and mass, then this amplitude allows us to separate the particles.
1 HDF, or the Hierarchical Data Format, is a commonly-used, cross-platform format for storing large amounts of data. C,

C++, Fortran, and Python libraries for reading it are available at https://support.hdfgroup.org/HDF5/.

103

All of the contaminants identified above have a larger nuclear charge Z than 46 Ar. Therefore, a good
first step to eliminating the contaminants is to remove any event that saturates the ion chamber ADC.
This was done by cutting any event where the digitized signal reaches zero. Next, events where the trigger
signal was not correctly digitized were removed, and so were events that did not trigger at least one CoBo.
This last criterion was found using a coincidence register recorded by the VME electronics.
After this first pass of cuts, the peaks were located in each ion chamber event and the peak heights
were recorded. Events with more than one peak were flagged as pileup and discarded. The peak heights
were then plotted in a histogram, shown in Fig. 5.8, and a cut was applied to keep only the large

46

Ar

peak. This peak contained 1.015 × 106 events. Note that the runs were divided into two sets in this plot
since the gain of the ion chamber was reduced between runs 216 and 217 by inserting a 0.3x attenuator
before the ADC. This was done to reduce the number of events that saturated the ADC.

5.3.4 Noise removal
In the next phase of the analysis, the merged events were read in from the HDF5 files produced in the
merging step. The baseline of each channel was corrected using the Fourier transform as described in
Section 3.1. Next, a single peak was found in each channel by finding the maximum, and the peak location was then refined by finding its center of gravity (see Section 3.2).
After peak finding, the (x, y, z, a) data points were calibrated as discussed in Section 3.2. This particular calibration used a nominal set of parameters and was required mainly to transform the z dimension
into spatial units; the data were re-calibrated later, so the specific values of these parameters were not
essential during this step of the analysis.
The calibrated data were processed using the Hough transform and nearest neighbor algorithms
from Section 3.3. This produced two values for each point: the minimum orthogonal distance to the nearest Hough line and the number of neighbors within a given maximum radius, which was set to 15 mm.
These two figures were then written to a new HDF5 file along with the uncalibrated peak locations in
three dimensions, the peak amplitudes, and the index of the pad that generated each point. Writing out
all of the data in this step allowed us to adjust our cuts interactively later, and to experiment with making
the cleaning cuts more or less restrictive as needed. In addition to this data, the center of curvature as
found by the Hough transform was also recorded in the output file.

104

5.3.5 Beam tracks
As discussed previously in Section 3.1, the beam tracks recorded during the

46

Ar experiment were not

usable for the analysis. This was caused in part by an unanticipated saturation effect in the preamplifier
of the GET electronics channels corresponding to the pads illuminated by the beam particle track. This
generated the unusually shaped signals shown in the top panel of Fig. 5.9. In addition, this saturation
also caused the pole-zero correction circuit to malfunction, causing the baseline restoration to fail. This
is shown in the inset in the top panel of Fig. 5.9.
Even if the electronics had worked as expected, the beam tracks would have still been difficult to
use in this experiment because of the size of their projection onto the pad plane. Due to an unfortunate coincidence, the drift velocity, electric and magnetic fields, and tilt angle used in this measurement
combined in such a way that the beam was projected onto a very small, almost point-like region on the
pad plane. The beam therefore illuminated very few pads, so the reconstructed beam tracks had very few
points, as seen in the bottom panel of Fig. 5.9. These sparse tracks are impossible to fit with a high degree
of precision.
These difficulties led us to abandon fitting the beam tracks entirely. Instead, we assumed that the
beam must lie along the z axis and introduced the vertex component of the objective function discussed
in Section 3.4.3.

5.3.6 Track fitting
The track fitting process began by reading in the cleaned data and center of curvature produced in the
previous step. A cut was applied to the data to keep only points that
• were less than 40 mm from the nearest Hough line,
• had more than two neighbors within a radius of 15 mm, and
• had a time bucket index (uncalibrated z position) less than 500.
The first two criteria were designed to exclude noise points as defined by the cleaning process. The final
criterion eliminated a narrow noise pulse that occurred in several channels at the end of many events,
perhaps due to switching noise in the electronics. Any event that had fewer than 50 data points remaining
after these cuts was discarded to avoid wasting computer time on extremely short tracks that correspond
105

Parameter

Value

Units

Tilt angle
Drift velocity
Electric field
Magnetic field
Gas pressure
CoBo write clock
Shaping time
Micromegas gain

6.2
ˆ − 0.552ˆv − 5.14w
ˆ
−0.0604u
ˆ
−1026ˆv + 9444w
ˆ
1.68w
19.2
12.5
280
500

deg
cm/µs
V/m
T
torr
MHz
ns
(unitless)

Table 5.1: Parameter values used while analyzing the 46 Ar(p, p ) data. The vectorial values are given in the
beam coordinate system of Fig. 3.2.
to very forward scattering angles in the center-of-mass system where the nuclear potential scattering is
minimal.
The remaining points were then calibrated and transformed into the beam coordinate system using
the method described in Section 3.2 with the parameters listed in Table 5.1. Initial estimates for the track
parameters were found using the method described in Section 3.4.4, and these values were used to seed
the Monte Carlo fitter. The experimental hit pattern was reconstructed using the peak amplitudes and
pad numbers from the calibrated data points. The optimizer was run on each track, and the results were
written to a set of SQLite2 databases for later analysis.
After track fitting was completed, these databases were collected into one dataset containing a total of 7.45 × 105 events. This is somewhat less than the number of

46

Ar events identified from the ion

chamber data since some very low-multiplicity tracks were eliminated after cleaning.
Poorly fit events were filtered out using a cut on the values of the objective function components
χ2 , as shown in Fig. 5.10. After the χ2 cuts, 2.2 × 105 events remained (approximately 30 % of the original
dataset). Further cuts were applied to remove events that were fit with a backwards scattering angle in the
laboratory frame (which is kinematically forbidden for elastic scattering), events that were not recorded
during a production run, and events that were marked as invalid due to misalignment between the GET
and VME DAQ systems (see Section 5.3.1). These cuts removed an additional 3 × 104 events, leaving a
total of 1.9 × 105 good events after cuts.
An interesting feature of both χ2 distributions is the very broad secondary peaks that occur at higher
2 A single-file relational database format. This provides the benefits of database transactions (such as preventing corruption

if the program is killed in the middle of a write operation) without the overhead and configuration requirements of a database
server.

106

values of χ2 . The events that make up these broad peaks consist partly of poorly fit tracks and noise
events caused by spurious triggers, as expected, but this region also contained many events from the
46

Ar beam particle scattering off of the carbon nuclei that make up part of the isobutane molecule. These

events were fit poorly by the optimizer since the algorithm assumed it was fitting the track from a proton
with Z /A = 1, and not, say, a 12 C nucleus with Z /A = 1/2. This implies that the χ2 cut can be used to filter
out carbon scattering events, a conclusion that was supported by extensive spot checks of the data.

5.3.7 Proton parameter distributions
Figure 5.11 shows the distribution of proton energies as reconstructed by the Monte Carlo optimizer.
The number of counts peaks at approximately 2 MeV and then diminishes rapidly at larger energies,
dropping below 10 around 6 MeV. The low number of counts at higher proton energies is explained by
the low trigger acceptance for these events, as shown in Figs. 4.7 to 4.11 and discussed in Section 4.2.
The distribution of reconstructed scattering angles in the center-of-mass frame is shown in Fig. 5.12.
The distribution is quite narrow, showing that the angular coverage is limited to angles between approximately 30° and 70° in the center-of-mass frame or 55° and 75° in the laboratory frame. Again, this is the
result of the low trigger acceptance outside this angular domain.
Figure 5.13 shows the distribution of vertex locations in the detector. The location (x 0 , y 0 ) of the vertex
was approximately normally distributed about the beam axis with a full width of approximately 10 mm
at half-maximum. Figure 5.13b shows that the distribution grows wider for smaller z 0 , which suggests
that this spread is caused by the emittance of the beam.

5.3.8 Vertex energy reconstructions
As discussed in Section 3.7, the 46 Ar vertex energy can be reconstructed for each event by two different
methods. One method uses the vertex location and a calculation of the beam particle’s energy loss in
the gas to calculate the final energy of the beam particle at the vertex location. The other method uses
the scattering angle and the energy of the outgoing proton to calculate the energy of the incoming 46 Ar
nucleus using formulas from relativistic kinematics. These two reconstructions are compared in Fig. 5.14.
It is apparent from this figure that the resolutions of the two reconstructions are vastly different: the
energy is reconstructed much more precisely from the vertex position than it is from the kinematics.

107

Runs 54-216

104

Count

103
102
101
100
10

1

Runs 217-269
103

Count

102
101
100
10

1

0

1024

2048

3072

4096

5120

Peak height [ADC bins]

6144

7168

8192

Figure 5.8: Distribution of peak heights in the ion chamber. The center of the 46 Ar peak is indicated by
the dashed vertical line, and the boundaries of the cut applied around this peak are shown with solid
lines.

108

Signal ampl.

3500
3000
2500
2000
1500
1000
500
0

100
50
0
0

256

Beam
Not beam

250
200

Uncal. y [mm]

512

150
100
50
0
50
0

100

200

300

ADC time bucket

400

500

Figure 5.9: A recorded beam track from the 46 Ar experiment. The top panel shows the traces from the lowgain region of the sensor plane, which is where the beam track is projected onto the pads. The signals
have an unusual shape which is much different from the signals generated by the proton tracks. The
inset illustrates the failure of the baseline restoration for one track. In the lower panel, the red squares
correspond to the reconstructed beam track, which is quite sparse.

109

12500

Count

10000
7500
5000
2500
0

0

20

40

2
pos

60

80

100

0

10

20

2
en

30

40

50

Figure 5.10: Distributions of the objective function components from fitting experimental data. The position component was cut at 14, while the energy component was cut at 10. Events with larger values for
either of these components were discarded.

104

Count

103
102
101
100
0

2

4

6

8

Proton initial energy (lab) [MeV]

10

12

Figure 5.11: Distribution of reconstructed proton vertex energies from the experiment. Note that the
vertical axis uses a logarithmic scale.

110

6000

Count

5000
4000
3000
2000
1000
0

0

30

60

90

120

Scattering angle (CM) [deg]

150

180

Figure 5.12: Distribution of reconstructed scattering angles in the center-of-mass frame.

111

20
500

15

400

5

Counts

y0 [mm]

10
0

300

5

200

10

100

15
20

20

15

10

5

0

x0 [mm]

5

10

15

0

20

(a) Distribution about beam axis

30
102

20

Counts

x02 + y02 [mm]

25

15

101

10
5
0

0

200

400

600

z0 [mm]

800

1000

100

(b) Spread as a function of z

Figure 5.13: Distribution of vertex positions in the detector. Plot (a) shows the overall distribution of
vertex positions in the plane transverse to the beam. Plot (b) shows that the radius of this distribution
increases for smaller z, or for vertex positions further from the beam entrance window.

112

4000

From position

3500

Counts

3000
2500
2000
1500
1000
500
0
4000

From kinematics

3500

Counts

3000
2500
2000
1500
1000
500
0

0

1

2

3

4

5

46Ar vertex energy (lab) [MeV/u]

6

7

Figure 5.14: 46 Ar vertex energies as reconstructed from the vertex position (top) and from kinematics
(bottom). The reconstruction methods are detailed in Section 3.7.

113

7

102

6
5

Counts

Vertex energy from kinematics (lab) [MeV/u]

8

4

101

3
2
1
0

0

1

2

3

4

5

6

7

Vertex energy from vertex position (lab) [MeV/u]

8

100

Figure 5.15: A comparison of the two vertex energy reconstructions. The line y = x is plotted for reference.

114

While the position-based vertex energy reconstruction appears to be quite precise, its accuracy depends on the precision of several experimental parameters as well as the accuracy of the SRIM [57] calculation used to generate the stopping power of the gas. The accuracy of the energy reconstruction can
be assessed by comparing the proton energy to the 46 Ar vertex energy. Assuming elastic scattering, kinematics calculations suggest that for proton energy E p and vertex energy E Ar ,
4E Ar ≈

Ep
2
cos (θ

lab )

.

(5.2)

These values are compared in Fig. 5.16, where the 2D histogram in the background shows values of
E p / cos2 (θlab ) and the red, dashed curve shows values of 4E Ar . This plot indicates that the calibration
was slightly off at low energies, so the data were recalibrated using the linear transformation
E (z) = E beam − α(E beam − E (z)),

(5.3)

where E beam = 4.17 MeV/u is the incoming beam energy, and α is a scaling parameter that controls the
magnitude of the recalibration. The best value of this scaling parameter was found to be 1.09, and the
recalibrated energy corresponding to this parameter is shown in Fig. 5.16 in cyan. This recalibration is
effectively the same as increasing the stopping power of the gas by 9 %, which indicates that the measured
gas pressure may have been somewhat inaccurate.
The poor resolution of the kinematic vertex energy reconstruction is caused by the fact that the beam
track was not usable in the recorded data. Because the beam track could not be fit, the emittance of the
beam introduces an uncertainty of up to 1° or 2° in the scattering angle since the analysis assumes the
beam lies along the z axis. Kinematics calculations show that this can introduce uncertainties on the
order of 1 MeV/u in the vertex energy at higher laboratory scattering angles (see Fig. 5.17).
Figure 5.15 shows a two-dimensional histogram comparing the two vertex energy reconstructions.
The distribution generally follows the line y = x, indicating that although the resolutions of the two reconstructions are vastly different, they tend to agree with each other on average. This suggests that the
data were correctly calibrated and that the fitting method was correct.
One prominent feature of the positional vertex energy spectrum is the sharp cutoff at 4.1 MeV/u.
This corresponds to the position of the beam entrance window at z = 1 m in the detector, and it is not
the product of any artificial cut applied to the data. Therefore, this edge can be fit to produce an estimate
of the resolution of the vertex energy reconstruction. This was done by fitting the edge with a modified

115

20.0

120

15.0

100

12.5

80

10.0
7.5

60

5.0

40

2.5
0.0

lab

0

200

400

600

Vertex position z0 [mm]

800

Count

lab) [MeV]

17.5

Ep/cos2(

140

Original
Recalibrated

20

65°
1000

0

Figure 5.16: Positional vertex energy recalibration. The 2D histogram shows values of E p / cos2 (θlab ), and
the two curves show values of 4E Ar . The original calibration is shown with a dashed red curve, and the
corrected calibration is shown with a solid cyan curve. The error band on the corrected calibration shows
the interval 1.065 ≤ α ≤ 1.115, where the best value is in the center at α = 1.09. Events with θlab > 65° were
eliminated from this plot to reduce the spread at higher energies caused by cos2 (θlab ) approaching 0.

3.00

5

2.75
2.50

Proton energy [MeV]

dEAr/d [(MeV/u)/deg]

4

2.25

3

2.00
2

1.75
1.50

1

1.25
0
50

55

60

65

70

Scattering angle (lab) [deg]

75

80

1.00

Figure 5.17: Derivatives of kinematic lines at a few proton energies. The vertex energy changes very
rapidly as a function of scattering angle as the scattering angle approaches 90° in the laboratory frame.

116

800

Fit
Data

700
600

Counts

500
400
300
200
100
0
3.8

3.9

4.0

4.1

4.2

4.3

4.4

Vertex energy from vertex position (lab) [MeV/u]

4.5

Figure 5.18: Results of a fit to the edge of the positional vertex energy distribution. The fit was performed
with the function defined in Eq. (5.4), yielding parameters µ = 4.11 MeV/u, σ = 19.7 keV/u, A = 694, and
b = 1.98. The bins in the energy dimension were 10 keV/u wide, and the error bars shown are ± N .
Gaussian cumulative distribution function (CDF)
f (x; µ, σ, A, b) = A(1 − Φ(x; µ, σ)) + b

(5.4)

where A corresponds to an amplitude, b is a vertical offset, µ is the mean, σ is the standard deviation,
and Φ is the usual definition of the Gaussian CDF, which is given by
Φ(x; µ, σ) =

x

1
2πσ2

−∞

2

e −(x −µ)

/2σ2

dx .

(5.5)

The result of this fit is shown in Fig. 5.18. The width was found to be σ = 19.7 keV/u, giving a resolution
of 46.4 keV/u at FWHM.

5.3.9 Cross sections and excitation functions
After reconstructing the

46

Ar vertex energies, the data were binned as a function of positional vertex

energy and center-of-mass scattering angle. The bins were chosen to be 50 keV/u wide in energy and 5°
wide in scattering angle. The binned results are shown as a two-dimensional histogram in Fig. 5.19, and
unnormalized excitation functions are shown in Fig. 5.20.

117

102

3

Counts

46Ar vertex energy (lab) [MeV/u]

4

2

101

1
0

0

20

40

60

80

100

Scattering angle (CM) [deg]

120

100

Figure 5.19: Two-dimensional histogram showing measured counts binned as a function of 46 Ar vertex
energy in the laboratory frame and scattering angle in the center-of-mass frame.
Normalized cross sections were found from this binned data using the general formula for a cross
section,
dσ
Nreac
=
,
d Ω Nbeam n t 2π sin θ∆θ
where Nreac is the number of scattering events that occurred, Nbeam is the number of

(5.6)
46

Ar nuclei that

entered the detector, n t is the number of target nuclei per unit cross sectional area in the detector, θ is
the scattering angle, which was taken in the center-of-mass frame, and ∆θ is the size of the angular bins.
The number of

46

Ar nuclei Nbeam that entered the detector was calculated using the ion chamber

data and the scalers:
Nbeam = NIC f Ar .

(5.7)

Here, NIC is the number of particles that passed through the ion chamber, a value that was recorded in a
scaler during the experiment. The correction factor f Ar is the proportion of the beam particles that were
46

Ar and not a contaminant. This was found by dividing the number of 46 Ar peaks in the ion chamber

data by the total number of peaks. This calculation does not account for particles that may have passed
through the ion chamber without being transmitted into the AT-TPC, but measurements of the transmission made during the experiment showed that this can be neglected.
The target density n t was found from the properties of the gas:
nt = n d x =

ρN A N t
M
118

1
dE.
d E /d x

(5.8)

Counts

CM = 40

800
600
400
200
0

Counts

600
400

400
200

200

0

0

CM = 55

CM = 60

CM = 65

300

400

200

200

100

0

0

200
150
100
50
0

CM = 70

125
100
75
50
25
0

CM = 75

20

40
20

10

0

0
CM = 90

CM = 95

10
8
6
4
2
0

15
10
5
0
0

1

2

3

CM = 80

30

60

CM = 85

Counts

CM = 50

600

600

Counts

CM = 45

800

4

5

Vertex energy [MeV/u]

10
8
6
4
2
0
0

1

2

3

4

5

Vertex energy [MeV/u]

0

1

2

3

4

5

Vertex energy [MeV/u]

Figure 5.20: Unnormalized excitation functions as a function of 46 Ar vertex energy in the laboratory
frame, shown at a variety of center-of-mass scattering angles. The error bars shown are ± N , where
N is the number of counts in a given bin.

119

In this equation, n is the number of gas molecules per unit volume, ρ is the mass density of the gas, N A
is Avogadro’s number, N t is the number of target nuclei per gas molecule, M is the molar mass of the
gas, d E /d x is the energy lost by the beam particle to the gas per unit track length, and d E is taken to be
the size of the energy bins. The energy loss data was calculated using SRIM, and the remaining quantities
were looked up in standard references.
The last remaining piece of Eq. (5.6) to calculate is the number of reactions Nreac . This quantity represents the total number of reactions that occurred in the detector, a value that is generally greater than the
number that were observed due to the trigger acceptance and dead time. This was calculated by correcting the observed number of counts using estimates of the proportion of events that were not recorded:
Nreac =
Here,

sim

Nobs
sim live pu corrupt

.

(5.9)

is the simulated total efficiency described in Section 4.2, which accounts for the acceptances

of the trigger and the analysis process. The second factor,

live ,

accounts for the dead time of the de-

tector, and was calculated as the ratio between the number of live triggers and the number of free triggers as recorded by scalers in the VME DAQ system. The factor

pu

accounts for the rejection of pileup

events, and was calculated by dividing the number of single-particle events in the ion chamber by the
total number of valid events. Finally,

corrupt

accounts for events that were discarded due to unrecover-

able misalignment between the GET and VME DAQ systems. This was calculated as a ratio between the
number of aligned events in each run and the total number of events.
To reduce the amount of uncertainty introduced by dividing the observed counts by the values of

sim

derived from the simulation, the simulated efficiencies were fit using kernel ridge regression, a machine
learning technique. Murphy [37] describes ridge regression as a modification of ordinary least squares
(OLS) regression that adds to the objective function a penalty term λ w
parameter λ and the square of the

2

2
2

depending on a weighting

or Euclidean norm of the fit parameters w. He states that this helps

keep the values of the fit parameters small, reducing the likelihood of overfitting. Kernel ridge regression
modifies this procedure further by replacing the inner products that appear in the OLS solution with calls
to a kernel function κ(x, x ) ∈ R that is generally nonnegative and symmetric in its arguments, and can
be interpreted as a measure of the similarity between feature vectors (measurements of the independent
variable) x and x [37]. The kernel function used in this fit was the radial basis function (RBF)
κ(x, x ) = e −γ
120

x−x

2

.

(5.10)

If the parameter γ is chosen to be σ−2 , this function corresponds to a Gaussian distribution [43]. The
fit itself was performed using a Python implementation of kernel ridge regression from the scikit-learn
package [43].
Figure 5.21 shows the result of applying this fitting procedure to the simulated efficiencies from the
third parameter set defined in Table 4.3. The fit curve, shown in orange, smooths out the small-scale statistical fluctuations in the efficiencies while maintaining the large-scale shape of the data. This smoothing is justified on the assumption that the efficiencies of the trigger and analysis should vary relatively
slowly as functions of energy. Therefore, using a fitted version of the data reduces the statistical noise in
the result without requiring an impractically large number of simulated tracks.
Figure 5.22 shows the cross sections calculated for the experimental data using Eq. (5.6). The measured cross sections are approximately of the same order of magnitude as the R matrix simulation result,
although the curvature of the data points is somewhat different from that of the simulation. This suggests that the simulation of the acceptance might be inaccurate, although this is difficult to prove due to
the low number of counts in each bin, especially at higher θCM .

5.3.10 Alternative analysis method
The accuracy of the cross sections presented in Section 5.3.9 depends heavily on the calculated efficiency.
This efficiency is found using a Monte Carlo simulation, so its value is subject to random fluctuations.
In addition, simplifying assumptions made in the simulation may introduce inaccuracies in the output.
Therefore, it is worthwhile to consider an alternative method of comparing the data to the R matrix
simulation.
In theory, the excitation function of a resonant scattering reaction should consist of isolated Breit–
Wigner resonances superimposed on a smoothly varying background described by the optical model [22,
pg. 597]. Thus, we should expect to see any resonances in the data as deviations from a smooth baseline
given by the optical model excitation function folded with the efficiency. If this baseline is removed, then
the resonances should be the only remaining features. This argument can also be applied to the results
of the R matrix calculation, yielding a theoretical curve that can be compared to the experimental data.
Since the analytical form of the baseline for the experimental data is unknown, it was modeled using
a quadratic function. Ordinary least squares regression was used to fit this function to the total recorded

121

CM = 40°

Efficiency

0.8

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0.0

0.0

0.0

CM = 55°

0.6

0.4

0.4

0.2

0.2

0.2

0.0

0.0

0.0

0.4

CM = 75°

0.20

CM = 80°

0.15

0.15

0.2

0.10

0.10

0.1
0.0

0.05

0.05

0.00

0.00

CM = 85°

CM = 90°

0.100
0.075
0.050
0.025
0.000

CM = 95°

0.04

0.06
0.04

0.02

0.02
0.00
0

1

2

3

CM = 65°

0.6

0.6

0.3

Efficiency

CM = 60°

0.8

CM = 70°

Efficiency

CM = 50°

0.6

0.8

Efficiency

CM = 45°

4

5

46Ar vertex energy [MeV/u]

0.00
0

1

2

3

4

5

46Ar vertex energy [MeV/u]

Fit

0

1

2

3

4

5

46Ar vertex energy [MeV/u]

Total efficiency

Figure 5.21: Simulated total efficiencies for parameter set 3 (as defined in Table 4.3) with kernel ridge
regression fit. Similar fits were performed for the efficiencies generated for the other parameter sets.

122

(d /d )CM [mb/sr]

CM = 40°

104

CM = 45°

6 × 103
4 × 103
3 × 103
2 × 103

103
6 × 102

(d /d )CM [mb/sr]

CM = 55°

2 × 103

CM = 50°

2 × 103

103

103

3 × 103

CM = 60°

CM = 65°

103

103

103
6 × 102
4 × 102
3 × 102

102
102

(d /d )CM [mb/sr]

103

(d /d )CM [mb/sr]

CM = 70°

104

102
101

CM = 75°

104
103
102
101
100

102
101
CM = 90°

104

102

101
2.0 2.5 3.0 3.5 4.0

EArlab [MeV/u]

CM = 95°

103

102

101

105
104

103

102

CM = 80°

103

CM = 85°

103

104

2.0 2.5 3.0 3.5 4.0

EArlab [MeV/u]

DSigmaIV

2.0 2.5 3.0 3.5 4.0

EArlab [MeV/u]

Experiment

Figure 5.22: Excitation functions in the region of interest for 46 Ar(p, p ) calculated from experimental
results. The results of an R matrix calculation performed with DSIGMAIV is also shown for comparison.
The parameters used in this calculation are given in Section 5.3.10.

123

1600

Counts

1500
1400
1300
1200
Quad. fit
Experiment

1100

2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00
46Ar vertex energy (lab) [MeV/u]

Figure 5.23: Experimental counts summed over all scattering angles and fit with a quadratic function.
The fit function was y(x) = −263(12)x 2 + 1687(74)x − 1207(108), where the numbers in parentheses are
the standard errors of the fit coefficients. A 1σ error band is shown around the fit using shading.
counts integrated over all scattering angles with 20 keV/u bins in the energy dimension. The result of the
fit is shown in Fig. 5.23. Integrating over all angles greatly increases the number of counts in each bin, and
it should not greatly affect the shape of the resonances since we are effectively only integrating between
approximately 30° and 60° due to the severe limitations on our acceptance. Nonetheless, any effect that
the integration has on the results will be compensated for below.
A similar procedure was applied to the R matrix calculation results from DSIGMAIV. Let R E θ represent
the calculated cross section in the bin defined by energy E and scattering angle θ. Then, a weighted total
excitation function was defined as
R˜E =

θ

R E θ Nθ
,
maxE (R E θ )

(5.11)

where maxE (R E θ ) denotes the maximum value of the cross section in each angular slice, and Nθ is a
normalization factor calculated from the experimental data. If the experimental bins are denoted D E θ ,
then these normalization factors are given by
Nθ =

DEθ
.
Eθ DEθ
E

(5.12)

This renormalizes the amplitude of the excitation function for each scattering angle to be proportional
to the number of counts observed in that angular bin in the experiment. Therefore, when the calculated
excitation functions are summed, the proportion contributed by each angular bin will be the same for

124

1

10

2

Full calculation
Baseline

Weighted total excitation
function [unitless]

10

2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00

Proton vertex energy (lab) [MeV]

Figure 5.24: Weighted total excitation function defined in Eq. (5.11) calculated for both the full R matrix
simulation (labeled “Full calculation”) and a calculation with no resonances (labeled “Baseline”).
the experiment and the calculation. One downside to this approach is that the resulting total excitation
function is unitless and its scale is arbitrary, but this will be addressed below.
Instead of fitting the R matrix calculation to remove the baseline, a second calculation was performed
using the same input file, but with all resonances removed. This produced the curve labeled “baseline”
in Fig. 5.24. This calculated baseline serves the same purpose as the baseline function that was fit to the
data above.
To put the data and the R matrix calculation on the same scale, a relative difference was calculated
between each quantity and its respective baseline. The relative difference was defined as
∆(D, B ) =

D −B
B

(5.13)

for data point D and baseline point B . This scales the two datasets to be of the same magnitude, and it
removes the baseline from each. The R matrix calculation curve was then convolved with a Gaussian of
standard deviation 19.7 keV to account for the resolution of the AT-TPC and provide a fair comparison.
The width of this Gaussian corresponds to the resolution of the positional vertex energy reconstruction
found in Section 3.7. The convolution greatly reduces the resonance amplitudes, especially in the case of
the narrow 3/2− resonance corresponding to the ground state of 47 Ar. This is shown in Fig. 5.25.

125

(value baseline) / baseline

0.04
0.02
0.00
0.02
0.04
0.06
0.08
0.10

Full resolution
Convolved with Gaussian
2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00

Proton vertex energy (lab) [MeV]

Figure 5.25: Effect of convolving the R matrix calculation with a Gaussian. The resonances in the convolved version are wider and have significantly smaller amplitudes.
5.3.11 Results
Figure 5.26 shows the final comparison between the R matrix calculation and the data. Four resonances
are shown in the result, and the properties of these resonances are given in Table 5.2. The R matrix calculation shown in the figure was performed with the Koning and Delaroche [29] global optical potential
parameters given in Table 4.2 on Page 62.
Each resonance was calculated with a resonance mixing phase of 20° to account for the fact that the
observed cross section is averaged over several fine structure states [47]. The value of this phase was
chosen to be compatible with calculations done by Harney [20] and with the results of Scott et al. [47]
with 40 Ar.
For each resonance, the total resonance width Γ can be written as
Γ = Γp + Γ↓ ,

(5.14)

where Γp is the proton single-particle width and Γ↓ is the spreading width caused by splitting over fine
structure states. The value of Γ↓ was taken to be 10(10) keV based on the values compiled by Reiter and
Harney [45] for proton scattering reactions on nuclei of a similar mass to 47 Ar. The relative uncertainty
of the values of Γp produced by DSIGMAIV was estimated to be 10 %.
CM
As discussed in Section 1.1.8, the resonance energy E res
is related to the energy E x of the excited state

126

0.100

(value baseline) / baseline

0.075
0.050
0.025
0.000
0.025

3
2

0.050

1+
2

1+
2

1
2

0.075
0.100

DSigmaIV
2.00

2.25

2.50

2.75

Experiment
3.00

3.25

3.50

46Ar vertex energy (lab) [MeV/u]

3.75

4.00

Figure 5.26: Comparison of data to R matrix calculation.

CM
E res

CM
E Ar

Ex

Jπ

S

Γ

2680(20)

−127 +39
−36 ± 20

0 +25
−20 ± 28

3/2−

0.27 ± 0.03 +0.21
−0.13

15(10)

2990(20)
3280(20)
3650(20)

183 +36
−50 ± 20

473 ± 36 ± 20
843 +32
−58 ± 20

310 +20
−40 ± 28

600 ± 20 ± 28
970 +10
−50 ± 28

1/2

+

1/2

+

1/2

−

0.30 ± 0.07 +0.14
−0.08
0.09 ± 0.02 +0.05
−0.07

0.42 ± 0.05 ± 0.09

30(10)
18(10)
34(10)

Γp
4.3(4)
20(2)
8.0(8)
24(2)

Γ↓
10(10)
10(10)
10(10)
10(10)

Table 5.2: Properties of the resonances shown in Fig. 5.26. Energies and widths are given in keV/u.
CM
E Ar
is the energy of the corresponding level in the analogue nucleus calculated assuming ∆EC − S n =
2807(30) keV, from the literature values of those parameters; E x is the excitation energy calculated from
the resonance assuming ∆EC − S n = 2680(20) keV, the energy of the ground state resonance; S refers to
the spectroscopic factor; Γ is the total resonance width; Γp is the proton single-particle width; and Γ↓
is the spreading width arising from mixing with the T< states. For quantities with two uncertainties, the
first value is the systematic uncertainty and the second is statistical.

127

in the analogue nucleus through the relation
CM
E x = E res
− (∆EC − S n ) ,

(5.15)

where ∆EC is the Coulomb shift of the state and S n is the neutron separation energy of the analogue nuCM
cleus. The energies E Ar
in Table 5.2 were calculated using a Coulomb shift of 6.471 39 MeV from a recent

theoretical calculation by Wang et al. [55]. To account for the fact that the Coulomb shift can vary slightly
between levels, it was assigned a systematic uncertainty of 30 keV, which was the standard deviation of
the differences between the observed and theoretical excitation energies found in 40 Ar [47, Table 1]. The
value of the neutron separation energy was 3.664 73(158) MeV, which was found using the NNDC Q-value
calculator3 , which is based on a recent set of mass measurements collected by Wang et al. [54]. In addition to the uncertainties arising from these values, a systematic uncertainty in the energy recalibration
(see Section 5.3.8) was added for each level. These uncertainties were [−20 keV, +25 keV] for the ground
state resonance, [−40 keV, +20 keV] for the first 1/2+ resonance, [−20 keV, +20 keV] for the second 1/2+
resonance, and [−50 keV, +10 keV] for the 3/2− resonance. These values were found by varying the recalibration scaling parameter up and down by 2.5 %, a limit which corresponds to the error band shown in
Fig. 5.16.
Alternatively, the excitation energy can be found by assuming that the ground state resonance must
correspond to zero excitation energy in the analogue nucleus. This produces the excitation energies
listed as E x in Table 5.2.
During the analysis, we noticed that the spectroscopic factors in the DSIGMAIV output were systematically much lower than the values found in the (d, p) measurement. This can be explained by the differences in the way the R matrix theory generates a single-particle width. As discussed in Section 1.1.6,
the R matrix model depends on a somewhat arbitrary boundary called the channel radius to set the scale
of the nuclear potential. Models of transfer reactions, on the other hand, do not include this division of
space, and the nuclear shape there depends only on the optical model parameters used. This produces
more realistic single-particle widths, and therefore more realistic spectroscopic factors. Thus, we renormalized the spectroscopic factors generated by DSIGMAIV to make them compatible with the values from
the optical model. These renormalization factors were found from an analysis of the 40 Ar case4 , and they
are given in Table 5.3. The spectroscopic factors listed in Table 5.2 are the renormalized values.
3 https://www.nndc.bnl.gov/qcalc/index.jsp
4 Analysis performed by Lisa Carpenter at the NSCL.

128

Jπ

Renormalization factor

3/2−
1/2+
1/2−

1.34(6)
1.22(24)
1.33(9)

Table 5.3: Renormalization factors applied to spectroscopic factors from DSIGMAIV.
Jπ

Slope [keV−1 ]

3/2−
1/2+
1/2+
1/2−

1.6 × 10−5
2.2 × 10−5
−7.7 × 10−7
1.3 × 10−5

I (Γ↓ = 10 keV)

Rel. uncertainty [%]

−1.55 × 10−3
−4.51 × 10−3
−4.28 × 10−3
−6.86 × 10−3

10
4.9
0.18
1.9

Table 5.4: Calculated values of the resonance integral from Eq. (5.17) and the relative uncertainty δI /I it
contributes to the spectroscopic factor for a variation of 10 keV in Γ↓ .
The systematic uncertainty in the spectroscopic factors can be related to variations in the integral
of each resonance with respect to energy. If we assume a Breit–Wigner resonance shape, then the cross
section has the general form
Γp
dσ
∝
.
dΩ
(E − E res )2 + Γ2 /4

(5.16)

The integral of this is
I≡

∞
−∞

Γp
(E − E res

)2 + Γ2 /4

dE =

2πΓp
Γ

,

(5.17)

which is inversely proportional to the spectroscopic factor since [22]
S∝

Γ
.
Γp

(5.18)

Thus, by measuring the change in I when varying Γp and Γ↓ , we recover information about the uncertainty in S. This was applied by varying Γ↓ from 0 keV to 50 keV in 5 keV steps and calculating the integral
of the negative part of the resonance5 at each step. The results of this are summarized in Table 5.4. The
dependence on Γ↓ was found to be small, so this effect was neglected in the final error analysis.
Finally, a Monte Carlo simulation was performed with TRIM [57] to assess the effects of energy straggling, or random fluctuations in the energy lost to the gas, on the energy resolution. The simulation
was performed using isobutane gas with density 6.108 × 10−5 g/cm3 (19.2 torr at room temperature). The
thickness of the gas target was varied between 10 cm and 70 cm in 10 cm steps, which corresponded to
5 The positive part depends more strongly on the resonance mixing angle, so it was excluded.

129

[keV/u]
en.

TRIM
Fit

11
10
9
8
7
6
5
4
1.5

2.0

2.5

3.0

46Ar vertex energy [MeV/u]

3.5

4.0

Figure 5.27: Results of energy straggling calculation with TRIM [57]. The orange line is a linear fit to the
data with slope −3.12(35) keV/MeV and intercept 16.2(10) keV. The shaded orange band is a 1σ error
band.
Energy (lab) [keV/u]

Straggling (lab) [keV/u]

2738
3055
3351
3729

7.67
6.68
5.16
4.58

Table 5.5: Straggling at the resonance locations.
vertex energies between 1.67 MeV/u and 3.86 MeV/u, thereby covering the range of interest. The vertex
energies were calculated by taking the mean of the outgoing

46

Ar ion energies, and the straggling was

taken to be the standard deviation of these energies. The results are shown in Fig. 5.27 with a linear fit,
and the straggling values predicted at the four resonance energies are shown in Table 5.5. The largest
predicted straggling was 7.67 keV/u, which is much smaller than the energy resolution or the energy bin
size, so the straggling was negligible.

5.3.12 Statistical testing
An F -test was performed to establish the statistical significance of these resonances. This test was constructed to compare two models to the data: the R matrix model and the horizontal line y = 0, which will
be referred to as the null model. The F -test asks whether the null hypothesis—that both models fit the
data equally well—can be rejected in favor of the alternative hypothesis—that the R matrix model fits

130

Jπ

Energy range [MeV/u]

F

p

3/2−
1/2+
1/2+
1/2−
—

(2.55, 2.86)
(2.86, 3.25)
(3.20, 3.55)
(3.4, 4.0)
(2.5, 4)

0.80
1.46
0.27
2.36
6.15

0.59
0.27
0.94
0.07
<0.01

Table 5.6: Statistical significance of the proposed resonances. The first three lines show the results of the
F -test applied to the three individual resonances. The last line shows the result of the test applied to the
data as a whole.
the data significantly better than the null model. This was done by computing the F statistic given by
F=

(RSS0 − RSSR )/(ν0 − νR )
.
RSSR /νR

(5.19)

Here, the subscript R refers to the R matrix model, and the subscript 0 refers to the null model. ν is the
number of degrees of freedom in each model, which was taken to be N − 6 for the R matrix model and N
for the null model, with N being the number of data points. Finally, RSS is the residual sum of squares
RSS =

i

ˆ i)
y i − y(x

2

(5.20)

ˆ
for data points (x i , y i ) and model prediction y(x).
This statistic was compared to the F distribution,
which has a cumulative distribution function given by
F (x; ν1 , ν2 ) = I

ν1 x
ν1 x+ν2

ν1 ν2
,
2 2

(5.21)

where I represents the regularized incomplete beta function. The null hypothesis was rejected if the pvalue
p = 1 − F (F ; ν0 − νR , νR )

(5.22)

was less than 0.10. This would imply that the R matrix model described the data significantly better than
a horizontal line, which would suggest that the resonance is statistically significant. Values of F and p for
each resonance are given in Table 5.6.

5.4 Discussion
5.4.1 Comparison to previous measurements
Table 5.7 and Fig. 5.28 compare the excitation energies E x and spectroscopic factors S found in this work
to literature values. In addition to the experimental values, three shell model calculations are included
131

in this table and figure. First is a calculation using the SDPF [41] interaction performed by Gaudefroy et
al. [18]. The SDPF interaction models the nucleus using a valence space consisting of the sd shell levels
for the protons and the p f shell levels for the neutrons [41]. This model is the oldest of the three, but
it provides the best agreement with our value of the energy of the 1/2− state. In addition to the SDPF
calculation, two calculations by Gade et al. [15] are included using the newer SDPF-U [40] and SDPF-MU
[52] interactions. According to Gade et al. [15], the SDPF-U interaction uses the same valence space as
the SDPF interaction, but it is tuned using experimental and theoretical data from neutron-rich silicon
isotopes to produce more-reliable results for these species. The SDPF-MU interaction improves this further by adding cross-shell tensor force interactions between the protons in the sd shell and the neutrons
in the p f shell, which helps the interaction to better reproduce B (E 2; 0+ → 2+ ) values and E (2+
1 ) and
E (4+
1 ) energies in silicon and sulfur isotopes near N = 28 [15]. These models agree quite well with the
experimental results of Gaudefroy et al. [16], but this is to be expected since the models were developed
using those results as input [15]. In general, the results of our experiment agree with the literature values
within 2σ for the ground state, but the agreement is less clear for the first excited state.
One possible explanation for the incompatibility between the 1/2− spectroscopic factor found in this
experiment and the larger value reported by Gaudefroy et al. [16] is the presence of a nearby 5/2− state.
This state was first observed by Bhattacharyya et al. [6] at 1234(4) keV, only 34(7) keV above the energy
of 1200(6) keV they observed for the 1/2− state, and it was later seen by Gade et al. [15] at 1231(4) keV.
This state was not known to Gaudefroy et al. [16], and with their experimental resolution of 75 keV for the
1/2− state, they would not have been able to separate this 5/2− state from the 1/2− state. Thus, we should
expect that their spectroscopic factor for the 1/2− state will be larger than the true value since their 1/2−
peak also includes contributions from transfer into the 5/2− state. The 5/2− state should not contribute
to the value of the spectroscopic factor of the 1/2− state observed in our experiment, however, since the
probability of populating this 5/2− state in a (p, p ) reaction is negligible compared to the probability in
a (d, p) reaction.
An interesting consequence of a smaller spectroscopic factor for the 1/2− state is its possible effect on
the single-particle energies. As discussed in Section 1.1.2, the single-particle energies are calculated from
a weighted sum of the experimental energy levels where the weights are the spectroscopic factors. Thus,
if the spectroscopic factor for the 1/2− state is substantially smaller than the value found by Gaudefroy et
al. [16], then the 2p 1/2 single-particle energy must be larger than the value they reported since it will have
132

3/2

3/2

20

0

Ex [keV]

20

40

1/2

900

0.00 0.25 0.50 0.75 1.00

S

1/2

1000

1100

Ex [keV]

1200

0.00 0.25 0.50 0.75 1.00

S

Present work
Gaudefroy et al.
Bhattacharyya et al.

SDPF-MU [Gade et al.]
SDPF-U [Gade et al.]
SDPF [Gaudefroy et al.]

Figure 5.28: Comparison between experimental results and literature values. Filled shapes indicate experimentally measured values, while empty shapes indicate shell model calculations.

Source
Gaudefroy et al. [16]
Bhattacharyya et al. [6]
SPDF [18]
SPDF-U [15]
SPDF-MU [15]

E x (3/2− )

S(3/2− )

E x (1/2− )

S(1/2− )

0
0
0
0
0

0.61(5)
—
0.64
0.711
0.634

1130(75)
1200(6)
1251
1139
931

0.81(6)
—
0.81
0.834
0.861

Table 5.7: Literature values of the measured excitation energies and spectroscopic factors. The SDPF,
SDPF-U, and SDPF-MU values are from shell model calculations performed by the authors of the referenced papers.

133

a smaller contribution from this lowest-lying 1/2− state and, therefore, a larger relative contribution from
higher-energy 1/2− states. If this single-particle energy rises without a change in the 2p 3/2 single-particle
energy, then the spin-orbit splitting of the 2p states will increase.
CM
Finally, the excitation energies reconstructed using the literature values of ∆EC and S n , shown as E Ar

in Table 5.2, are systematically lower than expected. The values in the column E x , calculated by taking
the resonance energy corresponding to the ground state as the value of ∆EC − S n , appear to be much
more realistic. This may indicate that the Coulomb shift calculated by Wang et al. [55] was larger than
the true value, which has not yet been measured experimentally.

5.4.2 Statistical significance
Of the four resonances described above, the only one that is statistically significant at the 10 % level is the
1/2− resonance corresponding to the first excited state of the 47 Ar nucleus. This resonance was the most
likely to be observed since, as a J = 1/2 state, it is expected to be broad compared to the higher-J resonances. The identification of the structure around 2.680 MeV/u as the ground state resonance is more
tenuous, however, since this structure was not found to be statistically significant. Finally, there appear
CM
to be two 1/2+ resonances at E res
= 2.990 MeV/u and 3.280 MeV/u, which correspond to center-of-mass

frame excitation energies of 0.310 MeV/u and 0.600 MeV/u in

47

Ar. These do not match any currently

known levels in 47 Ar. That being said, there is a 1/2+ state in 41 Ar at E CM = 1.88 MeV/u [47], and one was
also found in 45 Ar at 1.750 MeV/u [33], although these energies are substantially higher than the energies
reported above. Thus, it may be possible that these resonances correspond to previously unobserved
1/2+ states in 47 Ar that were not accessible in the previous (d, p) measurements. These states would also
be difficult to generate in a shell model calculation as it would require a valence space containing the
full sd and p f shells for both protons and neutrons, which would be very computationally intensive.
Alternatively, these resonances could also correspond to T< states in the compound nucleus, in which
case the true spectroscopic factors should be reduced by a factor of 2T0 +1 = 11 with respect to the values
reported in Table 5.2, which assumed a T> state.
A natural question about these results would be why the ground state resonance was not clearly
observed when the first excited state resonance was observed at a statistically significant level. There are
a few possible explanations for this. First, the ground state resonance is expected to be quite narrow as it

134

corresponds to a 3/2− state. The DSIGMAIV calculation predicted that this resonance should have a total
width of 15(10) keV, which is smaller than the energy binning used above and much smaller than the
FWHM energy resolution of the detector, which was found to be 46.4 keV/u in Section 5.3.8. This implies
that the resonance would be visible in at most 1-2 bins. Second, previous measurements using singleneutron transfer reactions [15, 16] have found that the spectroscopic factor for transfer into the ground
state is approximately 20 % smaller than for transfer into the first excited state. This implies that the
ground state resonance should have a smaller amplitude than the resonance corresponding to the first
excited state, which is already small. That being said, the results described above suggest that with better
statistics and the improved resolution that would be possible if the beam track could be reconstructed,
the ground state resonance should be observable with the AT-TPC.

135

CHAPTER 6
CONCLUSION

6.1 Summary
Studies of nuclei near the N = 28 shell closure can provide new insights into the structure of exotic nuclei
and the evolution of the shell model far from stability. To this end, an experiment was performed at the
National Superconducting Cyclotron Laboratory (NSCL) to measure resonant proton scattering on 46 Ar,
an N = 28 nucleus. This measurement was performed using the Active-Target Time Projection Chamber
(AT-TPC), a new detector built at the NSCL for experiments with low-energy, low-intensity radioactive
beams.
Since the AT-TPC outputs uncalibrated three-dimensional particle tracks for each recorded reaction,
a track reconstruction and analysis method had to be devised for this data. As described in Chapter 3,
the data were calibrated with simple linear transformations and then processed with an algorithm based
on the Hough transform to remove noise. Finally, a Monte Carlo algorithm was used to fit the tracks with
a model that includes contributions from the Lorentz force and the loss of energy to the gas that fills the
detector’s active volume. The location of each reaction vertex was determined by this fit, and the vertex
locations were combined with energy loss information to reconstruct the

46

Ar vertex energy for each

scattering event. This energy and the scattering angle determined by the fit were then used to produce
excitation functions.
After normalizing the excitation functions to units of cross section, they were compared to a simulation done with the R matrix code DSIGMAIV. The cross sections were of approximately the same order of
magnitude as the data, but the shapes of the distributions did not match well. We determined that this
was due to shortcomings in the simulation of the detector’s acceptance, so another analysis method was
devised (Section 5.3.10). This analysis revealed four possible resonances in the excitation function, one of
which was statistically significant. The significant resonance corresponded to the first excited 1/2− state
in 47 Ar, and the other three resonances corresponded to the ground state of 47 Ar and two 1/2+ states that
are either previously unobserved levels in 47 Ar or T< states in the compound nucleus.

136

6.2 Future outlook
As would be expected in a commissioning experiment, this measurement exposed a large and varied
collection of hardware problems, software bugs, analysis complications, and oversights in experimental
design that will need to be addressed in the future. Fortunately, many of these problems have already
been fixed or are in the process of being addressed.
Several of the largest problems we encountered in this experiment were caused by the large dynamic
range needed to resolve tracks left by both the heavy, highly ionizing

46

Ar beam particle and the light

scattered protons with Z = 1. To see the proton tracks, a large potential difference had to be applied
between the Micromegas mesh and the pads. This greatly over-amplified the beam tracks, causing an
unexpected saturation effect in the electronics that rendered these tracks unusable. This also created
the large baseline fluctuations described in Section 3.1. In future experiments, these problems will be
addressed with a programmable isolation circuit that will replace the simple circuit that is currently installed between the AsAd boards and the sensor plane. This programmable circuit will allow individual
pads to be biased to either a high or a low voltage, permitting the creation of a low-bias region for the
beam tracks. A prototype of this board was recently built and successfully tested, and a full complement
of these boards is now in the process of being procured.
Many of our remaining problems were caused by the trigger used during the experiment, which
greatly limited the angular and energy acceptance of the detector. This has been addressed in part by
creating a simulation of the trigger electronics (Section 4.1.3) that can be used before future experiments
to find appropriate trigger settings. In addition, some changes have recently been made to the CoBo
firmware to make the multiplicity trigger threshold setting more intuitive. These efforts, combined with
the availability of a true global multiplicity trigger from the now-available MuTAnT board, should make
triggering much more predictable in future experiments.
In addition to improving the trigger, having the full MuTAnT board fixes the DAQ system alignment
issues described in Sections 5.3.1 and 5.3.2. In future experiments, all CoBos will share a synchronized
global clock that can also be output to the VME DAQ system, allowing events to be merged much more
reliably by time stamp. This will also remove the ±1 time bucket uncertainty in z caused by the lack of
clock synchronization.
The last big problem we faced was the low experimental statistics. Somewhat ironically, one of the

137

largest limitations on the number of usable events was the high instantaneous beam rate injected into
the detector. This caused a substantial proportion (≥ 30 %) of events to be rejected due to pileup. This
could be addressed by accepting a lower beam rate, if possible. Additionally, the statistics for protonscattering experiments could be improved by filling the AT-TPC with pure hydrogen gas (H2 ) instead of
isobutane (C4 H10 ). This would eliminate the possibility of carbon scattering and increase the density of
proton scattering centers by a factor of 5 as compared to the same amount of isobutane gas. To make this
feasible even with the lower gain of pure hydrogen gas, a thick-GEM device has recently been simulated,
constructed, and installed in AT-TPC just in front of the Micromegas (see [46] for details). This provides
an additional degree of amplification that makes experiments in pure hydrogen and helium gas possible.
Notwithstanding the problems we encountered, as a commissioning experiment, this measurement
was largely a success. Although the resolution and statistics were not as good as anticipated, isobaric
analog resonances were visible in the recorded data, proving the feasibility of performing this type of
experiment with the AT-TPC in the future. Having addressed many of the limitations of this experiment,
future measurements should proceed much more smoothly and produce high-quality data.

138

BIBLIOGRAPHY

139

BIBLIOGRAPHY

[1] S. Anvar, P. Baron, B. Blank, J. Chavas, E. Delagnes, F. Druillole, P. Hellmuth, L. Nalpas, J. L. Pedroza,
J. Pibernat, E. Pollacco, A. Rebii, and N. Usher, “AGET, the GET front-end ASIC, for the readout of
the time projection chambers used in nuclear physic experiments”, IEEE Nucl. Sci. Conf. R., 745–
749 (2012) DOI:10.1109/NSSMIC.2011.6154095.
[2] M. Baranger, “A definition of the single-nucleon potential”, Nucl. Phys. A 149, 225–240 (1970)
DOI :10.1016/0375-9474(70)90692-5.
[3] P. Baron and E. Delagnes, AGET, a Front End ASIC for Active Time Projection Chamber: Data Sheet,
2013.
[4] S. Beceiro-Novo, T. Ahn, D. Bazin, and W. Mittig, “Active targets for the study of nuclei far from
stability”, Prog. Part. Nucl. Phys. 84, 124–165 (2015) DOI:10.1016/j.ppnp.2015.06.003.
[5] D. Bertsimas and J. Tsitsiklis, “Simulated annealing”, Stat. Sci. 8, 10–15 (1993)
1177011077.

DOI :10.1214/ss/

[6] S. Bhattacharyya, M. Rejmund, A. Navin, E. Caurier, F. Nowacki, A. Poves, R. Chapman, D. O’Donnell,
M. Gelin, A. Hodsdon, X. Liang, W. Mittig, G. Mukherjee, F. Rejmund, M. Rousseau, P. RousselChomaz, K. M. Spohr, and C. Theisen, “Structure of neutron-rich ar isotopes beyond N=28”, Phys.
Rev. Lett. 101, 032501 (2008) DOI:10.1103/PhysRevLett.101.032501.
[7] E. C. Bilpuch, in Isobaric spin in nuclear physics: proceedings of the conference on isobaric spin
in nuclear physics, edited by D. Robson and J. D. Fox (1966), p. 235.
[8] C. M. Bishop, Neural networks for pattern recognition (Clarendon Press, Oxford, 1995).
[9] B. A. Brown, Lecture notes in nuclear structure physics, Private communication, 2011.
ˇ
[10] V. Cerný,
“Thermodynamical approach to the traveling salesman problem: An efficient simulation
algorithm”, J. Optimiz. Theory App. 45, 41–51 (1985) DOI:10.1007/BF00940812.
[11] P. Descouvemont and D. Baye, “The R-matrix theory”, Rep. Prog. Phys. 73 (2010)
0034-4885/73/3/036301.

DOI :10.1088/

[12] Docker, https://www.docker.com/.
[13] R. O. Duda and P. E. Hart, “Use of the Hough transformation to detect lines and curves in pictures.”,
Commun. ACM 15 (1972) DOI:10.1145/361237.361242.
[14] R. Frühwirth, “Application of Kalman filtering to track and vertex fitting”, Nucl. Instrum. Meth. A
262, 444–450 (1987) DOI:10.1016/0168-9002(87)90887-4.
[15] A. Gade, J. A. Tostevin, V. Bader, T. Baugher, D. Bazin, J. S. Berryman, B. A. Brown, C. A. Diget, T.
Glasmacher, D. J. Hartley, E. Lunderberg, S. R. Stroberg, F. Recchia, A. Ratkiewicz, D. Weisshaar,
140

and K. Wimmer, “Single-particle structure at N=29: The structure of 47 Ar and first spectroscopy of
45
S”, Phys. Rev. C 93, 054315 (2016) DOI:10.1103/PhysRevC.93.054315.
[16] L. Gaudefroy, O. Sorlin, D. Beaumel, Y. Blumenfeld, Z. Dombrádi, S. Fortier, S. Franchoo, M. Gélin,
J. Gibelin, S. Grévy, F. Hammache, F. Ibrahim, K. W. Kemper, K. L. Kratz, S. M. Lukyanov, C. Monrozeau, L. Nalpas, F. Nowacki, A. N. Ostrowski, T. Otsuka, Y. E. Penionzhkevich, J. Piekarewicz, E. C.
Pollacco, P. Roussel-Chomaz, E. Rich, J. A. Scarpaci, M. G. St. Laurent, D. Sohler, M. Stanoiu, T.
Suzuki, E. Tryggestad, and D. Verney, “Reduction of the spin-orbit splittings at the N=28 shell closure”, Phys. Rev. Lett. 97, 092501 (2006) DOI:10.1103/PhysRevLett.97.092501.
[17] L. Gaudefroy, O. Sorlin, D. Beaumel, Y. Blumenfeld, Z. Dombrádi, S. Fortier, S. Franchoo, M. Gélin,
J. Gibelin, S. Grévy, F. Hammache, F. Ibrahim, K. W. Kemper, K. L. Kratz, S. M. Lukyanov, C. Monrozeau, L. Nalpas, F. Nowacki, A. N. Ostrowski, T. Otsuka, Y. E. Penionzhkevich, J. Piekarewicz, E. C.
Pollacco, P. Roussel-Chomaz, E. Rich, J. A. Scarpaci, M. G. St. Laurent, D. Sohler, M. Stanoiu, T.
Suzuki, E. Tryggestad, and D. Verney, “Gaudefroy et al. Reply”, Phys. Rev. Lett. 99, 099202 (2007)
DOI :10.1103/PhysRevLett.99.099202.
[18] L. Gaudefroy, O. Sorlin, D. Beaumel, Y. Blumenfeld, Z. Dombrádi, S. Fortier, S. Franchoo, M. Gélin, J.
Gibelin, S. Grévy, F. Hammaehe, F. Ibrahim, K. Kemper, K. L. Kratz, S. M. Lukyanov, C. Monrozeau,
L. Nalpas, F. Nowacki, A. N. Ostrowski, Y. E. Penionzhkevich, E. Pollacco, P. Roussel-Chomaz, E.
Rieh, J. A. Scarpaci, M. G. St Laurent, T. Rauscher, D. Sohler, M. Stanoiu, E. Tryggestad, and D.
Verney, “Study of the N = 28 shell closure in the Ar isotopic chain: A SPIRAL experiment for nuclear
astrophysics”, Eur. Phys. J. A 27, 309–314 (2006) DOI:10.1140/epja/i2006-08-047-0.
[19] Y. Giomataris, P. Rebourgeard, J. Robert, and G. Charpak, “MICROMEGAS: a high-granularity position sensitive gaseous detector for high particle flux environments”, Nucl. Instrum. Meth. A 376,
29–35 (1996) DOI:10.1016/0168-9002(96)00175-1.
[20] H. L. Harney, “Determination of spectroscopic factors from elastic proton scattering through the
isobaric analogue resonances”, Nucl. Phys. A 119, 591–608 (1968) DOI:10.1016/0375-9474(68)
90261-3.
[21] I. Heinze, “Development of a Hough Transformation Track Finder for Time Projection Chambers”,
Ph.D. (University of Hamburg, 2013).
[22] P. E. Hodgson, Nuclear reactions and nuclear structure (Oxford University Press, London, 1971).
[23] J. Illingworth and J. Kittler, “The Adaptive Hough Transform”, IEEE T. Pattern. Anal. PAMI-9, 690–
698 (1987) DOI:10.1109/TPAMI.1987.4767964.
[24] F. James and M. Roos, “Minuit: A System for Function Minimization and Analysis of the Parameter Errors and Correlations”, Comput. Phys. Commun. 10, 343–367 (1975) DOI:10.1016/00104655(75)90039-9.
[25] E. Jones, T. Oliphant, P. Peterson, et al., SciPy: Open source scientific tools for Python, http://www.
scipy.org.

141

[26] S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi, “Optimization by Simulated Annealing”, Science
220, 671–680 (1983) DOI:10.1126/science.220.4598.671.
[27] G. F. Knoll, Radiation detection and measurement, 4th ed. (John Wiley & Sons, Hoboken, NJ, 2010).
[28] T. G. Kolda, R. M. Lewis, and V. Torczon, “Optimization by direct search: New perspectives on some
classical and modern methods”, SIAM Rev. 45, 385–482 (2003) DOI:10.1137/S0036144502428893.
[29] A. J. Koning and J. P. Delaroche, “Local and global nucleon optical models from 1 keV to 200 MeV”,
Nucl. Phys. A 713, 231–310 (2003) DOI:10.1016/S0375-9474(02)01321-0.
[30] R. R. Labbe, Kalman and Bayesian Filters in Python (2015), https : / / github . com / rlabbe /
Kalman-and-Bayesian-Filters-in-Python/tree/360eccb.
[31] S. M. Lea, Mathematics for Physicists (Brooks/Cole – Thomson Learning, Belmont, CA, 2004).
[32] T. Lohse and W. Witzeling, “The Time Projection Chamber”, in Instrumentation in High Energy
Physics, edited by F. Sauli (World Scientific, Singapore, 1992) Chap. 2, pp. 82–156.
[33] F. Lu, J. Lee, M. B. Tsang, D. Bazin, D. Coupland, V. Henzl, D. Henzlova, M. Kilburn, W. G. Lynch,
A. M. Rogers, A. Sanetullaev, Z. Y. Sun, M. Youngs, R. J. Charity, L. G. Sobotka, M. Famiano, S. Hudan, M. Horoi, and Y. L. Ye, “Neutron-hole states in 45 Ar from 1 H(46 Ar, d)45 Ar reactions”, Phys. Rev.
C 88, 017604 (2013) DOI:10.1103/PhysRevC.88.017604.
[34] F. Marti, P. Miller, D. Poe, M. Steiner, J. Stetson, and X. Wu, “Commissioning of the Coupled Cyclotron Facility at NSCL”, in Cyclotrons and Their Applications 2001, Sixteenth International Conference (2001).
[35] K. I. M. McKinnon, “Convergence of the Nelder–Mead Simplex Method to a Nonstationary Point”,
SIAM J. Optimiz. 9, 148–158 (1998) DOI:10.1137/S1052623496303482.
[36] D. J. Morrissey, B. M. Sherrill, M. Steiner, A. Stolz, and I. Wiedenhoever, “Commissioning the A1900
projectile fragment separator”, Nucl. Instrum. Meth. B 204, 90–96 (2003) DOI:10.1016/S0168583X(02)01895-5.
[37] K. P. Murphy, Machine learning: a probabilistic perspective (MIT Press, Cambridge, MA, 2012).
[38] National Nuclear Data Center, NuDat 2 database, http://www.nndc.bnl.gov/nudat2/.
[39] J. A. Nelder and R. Mead, “A simplex method for function minimization”, Comput. J. 7, 308–313
(1965) DOI:10.1093/comjnl/7.4.308.
[40] F. Nowacki and A. Poves, “New effective interaction for 0 ω shell-model calculations in the sd − pf
valence space”, Phys. Rev. C 79, 014310 (2009) DOI:10.1103/PhysRevC.79.014310.
[41] S. Nummela, P. Baumann, E. Caurier, P. Dessagne, A. Jokinen, A. Knipper, G. Le Scornet, C. Miehé,
F. Nowacki, M. Oinonen, Z. Radivojevic, M. Ramdhane, G. Walter, and J. Äystö, “Spectroscopy of

142

34,35

Si by β decay: sd − f p shell gap and single-particle states”, Phys. Rev. C 63, 044316 (2001)
DOI :10.1103/PhysRevC.63.044316.

[42] D. R. Nygren, “The Time Projection Chamber — A New 4pi Detector for Charged Particles”, in
Proceedings of the 1974 pep summer study (1974), pp. 58–78, http : / / lss . fnal . gov / conf /
C740805/p58.pdf.
[43] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E.
Duchesnay, “Scikit-learn: machine learning in Python”, J. Mach. Learn. Res. 12, 2825–2830 (2011).
[44] E. Pollacco, S. Anvar, H. Baba, P. Baron, D. Bazin, C. Belkhiria, B. Blank, J. Chavas, P. Chomaz, E.
Delagnes, F. Druillole, P. Hellmuth, C. Huss, E. Galyaev, B. Lynch, W. Mittig, T. Murakami, L. Nalpas,
J.-L. Pedroza, R. Raabe, J. Pibernat, B. Raine, A. Rebii, A. Taketani, F. Saillant, D. Suzuki, N. Usher,
and G. Wittwer, “GET: A Generic Electronic System for TPCs for Nuclear Physics Experiments”,
Phys. Procedia 37, 1799–1804 (2012) DOI:10.1016/j.phpro.2012.02.506.
[45] J. Reiter and H. L. Harney, “Isospin mixing matrix elements extracted from isobaric analog resonances”, Z. Phys. A Atom. Nucl. 337, 121–129 (1990) DOI:10.1007/BF01294282.
[46] S. H. Rost, “Development of THGEM based detectors for AT-TPC applications”, M.S. (Michigan
State University, 2016).
[47] H. L. Scott, W. Galati, J. L. Weil, and M. T. McEllistrem, “Proton-Induced Reactions on 40 Ar and
Analogs of 41 Ar Levels”, Phys. Rev. 172, 1139–1148 (1968) DOI:10.1103/PhysRev.172.1139.
[48] A. Signoracci and B. A. Brown, “Comment on "reduction of the spin-orbit splittings at the N=28
shell closure"”, Phys. Rev. Lett. 99, 099201 (2007) DOI:10.1103/PhysRevLett.99.099201.
[49] O. Sorlin and M.-G. Porquet, “Evolution of the N = 28 shell closure: a test bench for nuclear forces”,
Physica Scripta T152, 014003 (2013) DOI:10.1088/0031-8949/2013/T152/014003.
[50] D. Suzuki, M. Ford, D. Bazin, W. Mittig, W. Lynch, T. Ahn, S. Aune, E. Galyaev, A. Fritsch, J. Gilbert,
F. Montes, A. Shore, J. Yurkon, J. Kolata, J. Browne, A. Howard, A. Roberts, and X. Tang, “Prototype
AT-TPC: Toward a new generation active target time projection chamber for radioactive beam experiments”, Nucl. Instrum. Meth. A 691, 39–54 (2012) DOI:10.1016/j.nima.2012.06.050.
[51] I. J. Thompson and F. M. Nunes, Nuclear reactions for astrophysics: Principles, calculation, and
applications of low-energy reactions (Cambridge University Press, 2009).
[52] Y. Utsuno, T. Otsuka, B. A. Brown, M. Honma, T. Mizusaki, and N. Shimizu, “Shape transitions in
exotic Si and S isotopes and tensor-force-driven Jahn-Teller effect”, Phys. Rev. C 86, 051301 (2012)
DOI :10.1103/PhysRevC.86.051301.
[53] A. C. C. Villari, D. Alt, G. Bollen, D. B. Crisp, M. Ikegami, S. W. Krause, A. Lapierre, S. M. Lidia,
D. J. Morrissey, S. Nash, R. J. Rencsok, R. Ringle, S. Schwarz, R. Shane, C. Sumithrarachchi, S. J.
Williams, and Q. Zhao, “Commissioning and First Accelerated Beams in the Reaccelerator (ReA3)

143

of the National Superconducting Cyclotron Laboratory, MSU”, in Proceedings of IPAC2016 (2016),
pp. 4–7, DOI:10.18429/JACoW-IPAC2016-TUPMR024.
[54] M. Wang, G. Audi, F. G. Kondev, W. Huang, S. Naimi, and X. Xu, “The AME2016 atomic mass evaluation (II). Tables, graphs and references”, Chinese Phys. C 41, 030003 (2017) DOI:10.1088/16741137/41/3/030003.
[55] Z. Wang, X. Guo, and R. Xu, “Effects of level inversion on the Coulomb displacement energies of
the Ar-K isobaric analog state pairs”, (2016), http://arxiv.org/abs/1605.08171.
[56] G. Welch and G. Bishop, An Introduction to the Kalman Filter, tech. rep. (University of North Carolina at Chapel Hill, Chapel Hill, NC, 2006), http://www.cs.unc.edu/~welch/media/pdf/
kalman_intro.pdf.
[57] J. F. Ziegler, M. D. Ziegler, and J. P. Biersack, “SRIM - The stopping and range of ions in matter
(2010)”, Nucl. Instrum. Meth. B. 268, 1818–1823 (2010) DOI:10.1016/j.nimb.2010.02.091.

144