Quick startup for Sunsort's Kinematics package group

This document gives a guide to starting to use Sunsort's Kinematics package quickly. It is aimed at users who are already familiar with using Sunsort to identify particles but want to churn through the kinematics quickly.

I'll assume you're analysing a standard CHARISSA experiment, in this case the reaction ¹²C(¹⁶O,¹²C,¹²C)⁴He.

Setting up the kinematics package for use

Getting access to the kinamtics package is easy. In FORTRAN simply add the lines

	include 'subs.i'
	include 'evsubs.i'
      

	#include <evsubs.h>
      

Telling Sunsort about the reaction

The first thing you need to do is to tell Sunsort about the reaction. The simplest way is to create a file called reac.cal in a directory called calib-inputs in the directory containing your sort code. The format if the file is documented in the documentation for the reac_init function call. So we do the following.

	% mkdir calib-inputs
	% cat > calib-inputs/reac.cal		(or use a text editor)
	65 1
	16 12
	12
	^D
      

We've chosen a dispersion of 1 that means we are working in MeV. We're not using the hitmass, but it has to be present.

Next, in your init code, add either the FORTRAN lines:

	if (reac_init(REAC_SORT).ne.SS_OK) then
	   print *,'Unable to initialise reaction'
	   call usererror(1)
	   return
	endif
      

or the C line

	if (reac_init(REAC_SORT) != OK)
	{
	    fprintf(stderr, "Unable to initialise reaction\n");
	    return 1;
	}
      

Initialising the event-struct array

You should next advise Sunsort how many particles you are dealing with so it can allocate sufficient memory. We'll tell Sunsort we want 7 particles (the 3 final state particles, the 3 possible reconstructed 2 fold particles and the 1 reconstructed compound nucleus). We'll tell it we want 8 detectors (not really important for this document) and we'll tell it the maximum reconstruction fold is 3. In FORTRAN this is done with

	if (ev_init(8,3,7).ne.SS_OK) then
	   print *,'Unable to initialise event-struct'
	   call usererror(1)
	   return
	endif
      

and in C by

	if (ev_init(evhead, reac, 8, 3, 7) != OK)
	{
	    fprintf(stderr, "Unable to initialise event-struct\n");
	    return 1;
	}
      

Per event initialisation

Each time you use the event-struct routines, it is important to clear out any data from the last event. Simply call ev_reset near the top of your sortin routine. In FORTRAN just do

	call ev_reset
      

and in C

	ev_reset();
      

Telling Sunsort about the particles

In your existing sort code, you are already identifying particles. Let's suppose the code has the following stages:

• For each detector segment (strip or quadrant)
  • See if anything happened in it, and if so
    • Identify the particle (its mass)
    • calculate its energy
    • calculate its direction of travel
    • Work out it's momentum by scaling its travel direction to sqrt(2*m*e).
• Rest of processing.

To tell Sunsort about the particles, we need at least the following information

• The particle's mass
• The particle's energy
• The particle's direction of travel

Which is conveniently exactly the information we have to work out with a standard sort code. Exactly how the particle's direction is specified depends on the sort code, but the Event-struct package will take it in a variety of formats. For the purpose of this example I'll assume you've calculated an unnormalised vector that points from the target to the position the particle hit the detector. I'll assume you have the particle energy in variable e, the mass in m and the vector in x, y and z.

In FORTRAN to tell Sunsort about a particle you do

	integer newp
	...
	newp=ev_new_particle()
	call evs_e(newp,e)
	call evs_m(newp,m)
	call evs_dir_xyz(newp,x,y,z)
      

in C it would be

	particle *newp;

	newp = ev_new_particle(evhead);
	evs_e(newp, e);
	evs_m(newp, m);
	evs_dir_xyz(newp, x, y, z);
      

These are the minimum parameters you need for the rest of the kinmatics, but you may want to add one or more of the following:

• A TDC signal in parameter t so you can later calculate time differences between particles.
• The particle charge in parameter z if you have Z resolution.
• The detector the particle was found in in parameter det.
• The detector segment (strip or quadrant) in paramter seg.

That was easy wasn't it.

Calculating event multiplicity

Once all the particles have beed found (that is, all the detectors have been looped over) you will want to check the event multipicity to make sure you detected the right number of particles.

In this example I'll assume we're going to detect the two ¹²C nuclei and not the ⁴He. So we want a something like this in FORTRAN

	if (evhead_n(1).eq.2) then
	   ... rest of analysis goes here ...
	endif
      

or in C

	if (evhead_n(evhead, 1) == 2)
	{
	    ... rest of analysis goes here ...
	}
      

Reconstructing missing particles

The next thing we need to do is to reconstruct our undetected particle by working out what momentum it was missing and creating a particle to carry that momentum. This is trivial with the Event-struct package. In FORTRAN we just do:

	j=missing_p()
      

where j is an integer variable. Any variable will do here we're going to throw the result away. In C we'd do

	missing_p(evhead);
      

and the result gets thrown away for us.

Etot calculation and Qggg gating

The next thing we need to do is calculate etot. This is also easy. Since we're going to use it a few times, we'll define a double precision variable called etot in the usual way, then we just do the following in FORTRAN

	etot = evhead_etot()
      

or in C

	etot = evhead_etot(evhead);
      

We will then want to plot etot, let's say into spectrum 100. In FORTRAN this is done with.

	call dinc1d(100, etot)		! or maybe etot*10
      

and in C with

	dinc1d(100, etot);		/* or maybe etot * 10 */
      

And then we want to gate on the Qggg peak for the rest of the analysis. Let's say we've put out etot limits into variables 1 and 2, then we do the following in FORTRAN

	if (etot.ge.var(1).and.etot.le.var(2)) then
	   ... rest of analysis goes here ...
	endif
      

or the following in C

	if (etot >= VAR(1) && etot <= VAR(2))
	{
	    ... rest of analysis goes here ...
	}
      

Excitation energies and resonant particles

This experiment can proceed via three ways (where C₁ and C₂ are the two ¹²C nuclei detected):

1. O + C -> C₁ + O^*-> C₁ + C₂ + He
2. O + C -> O^*+ C₂ -> C₁ + C₂ + He
3. O + C -> Mg^*+ He -> C₁ + C₂ + He

We want to calculate the excitation energies of the three possible intermediate (ejectile) nuclei. There are two ways to do this. We can calculate it purely from the recoil particle, or we can reconstruct the ejectile and calculate it from that. It turns out that in this reaction, the most accurate ways are to use the first method for the possible ¹⁶O^*ejectiles and to use the second method for the possible ²⁴Mg^*ejectile. I'll denote the three excitation energies in order as ex23, ex13 and ex23.

Since we know what particles we have and in what order, we can do this all very simply with (in FORTRAN)

	integer mg
	...
	ex23=-ev_q2(1)
	ex13=-ev_q2(2)
	mg=new_rp(1,2)
	ex23=ev_q(mg) + 13.93
      

or (in C)

	particle *c1, *c2, *mg;
	...
	c1 = evhead.ev + 1;
	c2 = evhead.ev + 2;
	ex23 = -ev_q2(c1);
	ex13 = -ev_q2(c2);
	mg = new_rp(c1, c2);
	ex23 = ev_q(mg) + 13.93;
      

We then have some complicated gating to do to determine which route we think the reaction took. Since this is a physics issue rather than a coding issue, I'll skip it.

Angular correlations

The final topic I'll cover here is that of angular correlations. For any resonant particle, you can calculate the angular correlation of its breakup. Let's suppose you want to make a theta*-psi plot in 2d spectrum 50 for the resonant ²⁴Mg particle shown above, this would be done in FORTRAN with

	call corr_polar(mg)
	call dinc2d(50,corr_psi(),corr_ths())	! maybe with some scaling
      

	struct correlation corr;
	...
	corr = corr_polar(mg);
	call dinc2d(50, corr_psi(corr), corr_ths(corr)) /* maybe scaled */