Exceeding the maximum beam scattering angle

[pvoss]

I have a question regarding an error output by Rachel when calculating simulated yields for the inverse kinematics Coulomb excitation of a 36Ar projectile on a 28Si target.

I have set up the calculation to integrate theta meshpoints (with endpoints defined by my particle detector's angular extent) between 10 and 30 degrees in the laboratory frame. From a kinematics spreadsheet calculator, I know that the maximum projectile scattering angle for this experiment at the “safe Coulex” energy is about 48.5 degrees in the laboratory frame (corresponding to a target recoil scattering angle of about 20 degrees). Thus, since in this inverse kinematics experiment the detector range spans the maximum scattering angle, I have input the target as the detected particle.

I've calculated the simulated yields for the logical experiments below (note that I've slowly decreased the minimum scattering angle to illustrate my soon-to-be-defined problem):
----------------------------------------------------------------------------
# 1 2 3 4 5 6
theta_min 15.0 21.0 17.0 15.0 13.0 11.0
theta_max 42.0 30.0 30.0 30.0 30.0 30.0
phi_min 0.0 0.0 0.0 0.0 0.0 0.0
phi_max 360.0 360.0 360.0 360.0 360.0 360.0
un-inv. Z 14 14 14 14 14 14
un-inv. A 28 28 28 28 28 28
Ebeam 63.0 63.0 63.0 63.0 63.0 63.0
Safe energy 64 64 64 64 64 64
Tgt. thickness [mg/cm^2] 0.6 0.6 0.6 0.6 0.6 0.6
Detected t t t t t t
Inverse kin. True True True True True True
Kinematic solution -- -- -- -- -- --
Norm to expt # 1 2 3 4 5 6
Somm. param. 33.0 33.0 33.0 33.0 33.0 33.0
# of Ge detectors 1 1 1 1 1 1
--------------------------------------------------------------------------------

I noticed that even though all of the theta_min and theta_max ranges for scattered target nuclei span the kinematics solution corresponding to the maximum projectile scattering angle (except Experiment 2), I only get the following error for Experiment 6:
--------------------------------------------------------------------------------
Gosia reported an error:
ERROR- MAXIMUM SCATTERING ANGLE IS 48.00 DEGREES FOR EXPERIMENT 6

********* END OF EXECUTION **********

---------------------------Error----------------------------

You have exceeded the maximum beam scattering angle
(calculated inelastically) for this experiment. There are
two ways that you might fix this problem: (1) Define the
experiment(s) so that the target particle is detected in the
case of inverse kinematics, or (2) Set the state for
calculation of the scattering kinematics to the ground state
using the Gosia Controls button, option k.

Done.

Fri Sep 30 13:11:01 PDT 2011 ...
Processing cannot continue properly due to the error above.
--------------------------------------------------------------------------------

I get this error for the same experiment (11deg<theta_lab<30deg) even if I shuffle their order, so it's not tied to the last experiment listed.

If I understand correctly, I have properly defined the experiments such that “the target particle is detected in the case of inverse kinematics” and therefore I shouldn't get this error? Or perhaps I am making a silly mistake. Either way, a hint would be helpful so I know whether I can put faith in my calculations or not.

Many thanks and cheers,
Phil

[hayes]

Hi Phil,

Can you post the level data, so that I have the right Q-value? I will try to reproduce the error and fix it. The easiest way for me is in the "rachel" format. If you don't have that already, you can go to Tools-->"el" (export level scheme), save to a new file.

Also, I will have a new version ready in about a week (or ask me if you want to preview it) that will address bugs and have a better simulation output with understandable efficiencies instead of the poorly-defined "intrinsic efficiency" in the beta version.

Adam

[pvoss]

Hello Adam,

Here is the level data:
---------------------------------------------------
18 36
gsb 0.0 + 0.0
gsb 2.0 + 1970.39
gsb 4.0 + 4414.36
-------------------------------------------------

I tried attaching the pickle.jar file too, but was not allowed to upload a file with ".jar" extension. Thank you very much for your help.

Cheers,
Phil

[hayes]

Hi Phil,

I found a simple solution. If this is not acceptable, or if you think there are other accuracy questions, I can see if there is another way to avoid the error.

The problem comes partly from the inelastic kinematics using the Q-value of the 2+ state. (By default, Gosia and Rachel use the second state in the level scheme to get estimate the Q-value. It is necessary to pick a single Q-value for each calculation in a semi-classical code, but the calculation can be repeated, of course, to use different Q-values for each state.)

I got the simulation to run without errors by changing from the default Q-value of level #2 to the ground state (Q=0).

This is done by clicking "Gosia controls" at the left and following the prompts as below:

Code: Select all

Menu:
Physical parameters
   v Vacuum deorientation control
   k State for calculation of the scattering kinematics
Data handling options
   t Normalization transition for all normalized yields **new**
   n Data set normalizations
   w Data set weights  **new**
Accuracy and calculation speed options
   i E-theta integration accuracy
   f Options for fitting of matrix elements
Enter choice (one character) from above: k
Current state for calculation of kinematics is  Band " gsb ", spin  2.0
Band NUMBER of new kinematics state: 1
Spin: 0
Kinematics state set to  Band " gsb ", spin  0.0


There have been a few questions raised recently about using Q=0 for the scattering cross section, and repeatedly calculating yields of other states using their own Q-values. You can find more about this on the forum.

At the moment, for your 5th experiment, there is about a 10% change in the 2+ --> 0+ yield between the two Q-values. Obviously, this must be addressed in analysis, but this solution may be accurate enough for your simulation.

Please let me know.

Adam

[warr]

I think the problem is occurring when you happen to have a meshpoint at a target angle of 20.5° (exactly halfway between 11 and 30°). With your kinematics the maximum projectile angle should be almost exactly 48°, which corresponds to a target angle of 20.5° (this is calculating for the exit energy. Obviously it is more like 48.5° for the incident energy, but from your error message, it must be the exit energy which is triggering the error).

The code takes the angle of the meshpoint (which is a target angle) and converts using INVKIN into projectile angles. Then CMLAB calculates the maximum projectile angle. In principle, this should never fail. In fact, if there is a small rounding error, this could lead to angles very close to the maximum being considered beyond the maximum.

I wonder if it is possible to tweak the limits in your input (i.e. a small change) so it doesn't calculate a meshpoint at exactly the critical angle. I'm not sure how much would be needed, but perhaps not enough to make any real difference to the calculation. In which case that would be a better workaround than changing NCM, which changes the kinematics (i.e. the kinematics is calculated based on the ground state not the first excited state). Also I think it is by chance that this works, because it changes the floating point values just enough not to have the rounding error problem.

So this would seem to be a bug in my code (INVKIN). I'll have to think about how to fix it. Probably it will involve trapping the case of being only a little over the maximum angle in CMLAB which is part of the original Gosia.

Thank you for point it out.

[hayes]

Hi Phil and Nigel,

Nigel's suggestion to tweak the meshpoints seems to work. Your (Phil's) input had 21 theta meshpoints for the problematic experiment (#6). I changed it to 18 meshpoints, and it worked. It also worked for 24 meshpoints and gave an identical output.

This is a better solution than being forced to change the Q-value. This way seems works for the Q of the 2+ state, and the integration is obviously converged, since I get identical output to 5 digits with 18 or 24 meshpoints.

I did this through the GUI. Here's how:
Click "Gosia controls", then follow these prompts.

Code: Select all

Menu:
Physical parameters
   v Vacuum deorientation control
   k State for calculation of the scattering kinematics
Data handling options
   t Normalization transition for all normalized yields **new**
   n Data set normalizations
   w Data set weights  **new**
Accuracy and calculation speed options
   i E-theta integration accuracy
   f Options for fitting of matrix elements
Enter choice (one character) from above: i
The accuracy can be increased at a small expense of speed
during the "Integrated yields" operations by increasing the
number of interpolation points for the beam energy, polar
scattering angle (theta) mesh.

The accuracy can be increased at a greater expense of speed
by increasing the number of energy and theta meshpoints.

NOTE: If you decrease the number of interpolation points,
or meshpoints, you should check that the integrated
(calculated) yield accuracy has not been reduced.  See
"integration accuracy" in the Help function.

Present E-theta meshpoint settings:
                 beam    scattering
                 energy  theta
-----------------------------------------------
Experiment  1 :    20      30
Experiment  2 :    20      20
Experiment  3 :    20      20
Experiment  4 :    20      20
Experiment  5 :    20      20
Experiment  6 :    20      21
Would you like Rachel to choose the optimum settings for all experiments [y/N]?
Experiment number to change, "all", or return to quit:6
Enter one integer for energy meshpoints and one for theta on the same line. 20 18
                 beam    scattering
                 energy  theta
-----------------------------------------------
Experiment  1 :    20      30
Experiment  2 :    20      20
Experiment  3 :    20      20
Experiment  4 :    20      20
Experiment  5 :    20      20
Experiment  6 :    20      18
Experiment number to change, "all", or return to quit:
Quitting.
Present settings of integration subdivisions for interpolation.
                 beam    scattering
                 energy  theta
-----------------------------------------------
Experiment  1 :    50      50
Experiment  2 :    50      50
Experiment  3 :    50      50
Experiment  4 :    50      50
Experiment  5 :    50      50
Experiment  6 :    50      50
Would you like to set the default number of meshpoints for all experiments [y/N]?
Experiment number to change, "all", or return to quit:
Quitting.


Adam

[warr]

I have found an error in my INVKIN subroutine. The variables tau and taup are not declared explicitly, so they are implicitly REAL*4. However, the equivalent variables in CMLAB are REAL*8. This means that the precision in INVKIN is less than in CMLAB. This won't make any noticeable difference to the result, but it can happen, as in your case, that it triggers the problem you have.

So the simple solution is to declare these two variables in the subroutine INVKIN as REAL*8.

Unfortunately, that is not the whole story. It is still possible to find cases where the difference between INVKIN and CMLAB (which are essentially a function and its inverse) still lead to small differences that can push the value of tau * SIN(theta_p/57.297795) above unity (we are talking of around 2 * 10^-16 above unity here). Since Gosia has to calculate the arcsine of this quantity, it will fail. This goes deep into the limitations of floating point arithmetic and is not easy to solve. So if other people hit this problem, the solution (after declaring tau and taup properly) is to tweak some parameter by a small amount. I found that changing the excitation energy by 0.01 eV (eV not keV!) was enough in the case I found. Fortunately, such cases are quite hard to find even when you are looking for them, so you have to be very unlucky to find them (after fixing the precision in INVKIN, that is).

[hayes]

Note: At the moment we don't have the rachel patches in the version control system that Nigel manages. You *can* however use any recent version of Gosia with rachel, as long as you don't try to use the functions in the GUI that plot the Coulex amplitude or collision function. These functions are for troubleshooting only. The calculations we are talking about can be repeated with a new version, but there should be no difference since the integration is obviously converged in this case. (I tested this by changing the number of meshpoints in theta up and down and finding no difference in the yields reported.)

Adam

[pvoss]

Adam and Nigel,
Thank you very much for your help!
Phil