Gosia2 brute force error estimation

[C.Bauer]

Hi to all,
I used the Gosia2 error method as described in the Wiki:
http://www-user.pas.rochester.edu/~gosia/mediawiki/index.php/Error_estimation

For some reason I am not so confident that I am doing everything correctly:
1. The new minimum of the parabola is ~0.815 for the investigated M.E. while it was ~0.835 before after the regular Gosia2 minimization, should I worry about that deviation?
2. The new error is smaller (~0.02) compared to ~0.06 which the standard error calculation from Gosia gives. But I thought the reason I am doing all this is because otherwise I UNDERestimate the error...

Some more information: It is MINIBALL data, 140Nd -> 48Ti, and I just see the 2+1 -> 0+1 transition in both isotopes. The scattering angle was divided in 2 "experiments".

The described procedure says in step 3.1 "Calculate the corrected yields using OP,CORR." So I use NTAP=3 and OP,INTI, correct? In step 3.2 I minimize again with NTAP=4 and OP,MINI. Should I put OP,REST also? Do I also have to use the OP,MAP command in the iterations? These are things I am not sure about...

Finally, the reduced chi^2 values from the output files I multiply by 2 (2 experiments, 1 yield in each) and add them up to get the chi^2_total, which is the one I have plotted.

That's it so far, I would be very thankful for any hint or any comment if you had a similar experience.

Best regards,
Christopher

[hayes]

Hi, Christopher.

I have a feeling we may need to mail back and forth a few times to sort this out. Would you like to send me your set of files? One thing I could do is to "simulate" data as I like to do and see what errors I get when I know the answer. You can do this yourself by running OP,INTI and putting the output into the yield files (tape 3) as though they were data. You should get a chi-squared of zero, but you can still find the errors using each method. Using fabricated data is a good way to test methods.

1. The new minimum of the parabola is ~0.815 for the investigated M.E. while it was ~0.835 before after the regular Gosia2 minimization, should I worry about that deviation?

The minimum will usually change by a small amount in any method of error estimation. I think this is probably because of the step size chosen by Gosia around the minimum. I wouldn't worry about this if it doesn't change by much more than the error, which seems to be the case here.

2. The new error is smaller (~0.02) compared to ~0.06 which the standard error calculation from Gosia gives. But I thought the reason I am doing all this is because otherwise I UNDERestimate the error...

Yes, I think the error reported by Gosia should be smaller, if there is any difference. At least that is my limited experience with Gosia 2, and it seems to be generally the case for Gosia 1. As Nigel told you, maybe it won't always be the case. At the moment I don't understand why you get a smaller error this way. Are you re-running the fit (OP,MINI) for both inputs at each step in the brute-force method? I have a feeling that you understand the procedure as well as I have it written out.

Have you tried changing the OP,MINI settings, like the chi-squared goal and the convergence criterion? If it reports that it is stopping the fit because it reached the chi-squared goal, you might try lowering it and see what happens if it terminates because of convergence.

Gosia's reported reduced-chi-squared should be for the data points (total yields for all detectors and experiments) in one input as you said. For example, if the reduced chi-squared in the output file for the target is 1.0 and there are 2 data points for that input, then the total chi-squared for that input is 2.0. This will be added to the total chi-squared for the projectile to get the total chi-squared. Is this the way you have been summing it?

The described procedure says in step 3.1 "Calculate the corrected yields using OP,CORR." So I use NTAP=3 and OP,INTI, correct? In step 3.2 I minimize again with NTAP=4 and OP,MINI. Should I put OP,REST also? Do I also have to use the OP,MAP command in the iterations? These are things I am not sure about...

This is correct. Just to be sure, in the OP,CORR step you are putting OP,CORR just after OP,INTI and before OP,EXIT each time you iterate the corrected yields, right? This seems obvious, but I just wanted to be sure. Then you change NTAP back to 4 for OP,MINI. Again, this is probably obvious at this point.

Let me know what you think so far, and I'll keep thinking about it.

Adam

[C.Bauer]

Dear Adam,
thanks for the quick answer. Yes, so I think I am doing the correct steps and also the correct summing of the chi-squares. We have a conference right now in Darmstadt, after that I will try to change the criterion and the other ideas you suggested. Nevertheless I can send you my files by mail, so you can have a look on your own and maybe there is an obvious mistake I was not aware of...
Christopher

[C.Bauer]

I have to add something:

I did the same "brute force" minimizations again, now keeping 2 matrix elements fixed and calculating chi^2. The 2 matrix elements were M_20 and M_22, that means the diagonal connected to the quadrupole moment of the 2+1 also. This gives me a chi^2 map like in this nice new publication http://prl.aps.org/pdf/PRL/v108/i6/e062701.

Now the error is much bigger, when projecting the chi^2+1 contour on the corresponding axis! In my opinion I should get something similar, which tells me there is still something not consistent behind what I am doing...

M_20 = 0.845 (+/- 0.055) roughly

[hayes]

Hi Christopher,

I hope this is more helpful than it is long!

The first numbered items are very important. After these, I have other
suggestions, but the first ones must be addressed first, or your results may
not be accurate.

1. I see that you have all experiments normalized to experiment 1 within each
input. It seems that experiments 2 and 3 in the 140Nd input are at a different
beam energy from experiment 1. If you try to normalize an experiment to one
from another beam run, you will get ridiculous results, unless you have some
very good charge integration. If you do, please let me know so that I
understand.

2. I see that you have target-detection experiments normalized to other
experiments. This will not work in Gosia, but there will be no complaint in
the output. (This is in the new manual, but we have not added an error in the
code.) If you normalize a target-detection experiment to any other target- or
projectile-detection experiment, you will get systematic errors in chi-squared
that can be huge in some cases.

3. If you are able to accurately change these experiments to the equivalent
projectile scattering range, then you can normalize them together. If you do,
read the newest Gosia manual (OP,YIEL section) to see how to calculate the
normalization constants. Even if you measure them, you will need extra factors
that are not intuitive (but easy to calculate).

4. Be careful if you change to the equivalent projectile scattering range. I
haven't looked at your kinematics, but this can present accuracy problems,
which are described in the scattering theory section of the new manual. If you
don't want to do this, then normalize each experiment individually by setting
LN to its own number:

EXPT
3,60,140
-30,64,335.,20.0,3,1,0,0,360,1,1
-30,64,366.,-26.0,3,1,0,0,360,1,2 ! was 1
-30,64,366.,-23.0,3,1,0,0,360,1,3 ! was 1

Since gosia2 normalizes beam to target excitation, you will still have enough
sensitivity to the matrix elements, if you don't have too many fit parameters.

----------------------------------------------------------------------------------

The comments below may help, but the ones above should be considered/fixed before trying anything below.

I had trouble running your inputs. I think these files are not all from the
same calculation, because they reference 140Nd_Zn* and the yield file has only
2 experiments, while EXPT lists 3. But that is ok. I now see that it would
take a lot of time to repeat what you're doing, so I'll comment on what I read
in the input files. I may not have too much time to re-run more than a few
calculations.

Fit parameters:

I see 6 total free parameters in the inputs, and 8 data points (5 experiments
plus 3 previous measurements). I think that even in gosia2 the overall
normalization could be an additional parameter, but I sometimes get confused
about this. If so, then you have 7 free parameters in total. I'll think about
this again later. Usually, you need many more than 7 data points to fit 7
parameters, but with the high precision in the previous lifetime measurements,
it seems like 6 user-parameters might be ok.

My confusion is that with too many fit parameters, it seems like you would get
a much *larger* correlated error, since the additional free parameters should
overcompensate for a change in one individual parameter.

Maybe the following questions will help to clear this up:

Do you see major conflicts between some of the data? For instance, do you see
that a lifetime and a fitted matrix element are in large conflict?

You could try a test fit with more parameters locked. Do you get a sensible
correlated error in the one-dimensional method, and roughly the same error in
the 2D method? I get a lot of insight by calculating yields with OP,INTI and
then putting them back in the *.yld files to do test fits. If you can't do a
fit with simulated data, then you almost certainly will have a problem with the
real data, using the same fit parameters, etc.

Can you take a few more points around the minimum in the 1-D error calculation?
I see only 1 point within the 1-sigma limits on the pdf plot, if I understand
correctly. The shape may have a smaller slope around the minimum, and going up
to chi-squared = 10 probably doesn't help you to understand the shape around
the minimum.

How are the fits terminating? Sometimes OP,MINI parameters need to be
adjusted, such as if the fit terminates because it hit chi-squared=1. and the
best fit value is a little less than 1.

-----------------------------------------------------------------------------------

Some other concerns:

It seems like you have data only for the 2-->0 transition, and you're using the
4+ state as a "buffer state." This is good.

You have a weight of 10 for the spectroscopic data (lifetimes, etc.). Do you
mean to weight them so heavily? This adjustable weight can be somewhat useful,
but I don't think that there is a real mathematical justification to pick any
number other than 1. The error bars should take care of the weighting, in my
opinion; changing the weight seems no different from changing the error bars,
and that should be done individually for each measurement, if at all. If there
is no justification I haven't thought of, then increasing their weight
artificially reduces the correlated error measurement.

You could try using OP,MAP every time you use OP,MINI. This will update for
the present set of matrix elements. I have found that this can be done in the
same step, preceding OP,MINI with OP,MAP:

OP,REST
0,0
OP,MAP
OP,MINI
2100....

It would be worth trying to run OP,CORR (with OP,INTI and NTAP = 3) at each
minimization step. Sometimes this is not too important, but it will be if the
correction factors change significantly as a function of the matrix elements.

Adam

[hayes]

Hi Christopher,

One last thing I noticed--you have 3 experiments in the projectile input and 2 in the target input. (Maybe this is because of a file mix-up.) You may not get an error in Gosia2, but I don't think you will get sensible results in a fit if you don't have the same experiments in the same order in both inputs.

Adam

[C.Bauer]

I am very sorry, my mistake! The code at some point has overwritten the input file from runs with a 2nd target, namely Zn-64. That's not affecting the results I reported before and in my email. For now I just want to do the analysis with the Ti-48 target, I send it to you for completeness, so I have 2 experiments defined in both, projectile and target, of course.

About the other issues I will report later...

[C.Bauer]

1. I see that you have all experiments normalized to experiment 1 within each
input. It seems that experiments 2 and 3 in the 140Nd input are at a different
beam energy from experiment 1. If you try to normalize an experiment to one
from another beam run, you will get ridiculous results, unless you have some
very good charge integration. If you do, please let me know so that I
understand.

The experiments are all from the same beam run with the same energy. I don't understand completely what to put as mean energy and mean angle. Of course, they have to be in the corresponding range. The reason why I used different energies is that once somebody told me to vary them a little bit to get correction coefficients close to 1 after the integration.

2. I see that you have target-detection experiments normalized to other
experiments. This will not work in Gosia, but there will be no complaint in
the output. (This is in the new manual, but we have not added an error in the
code.) If you normalize a target-detection experiment to any other target- or
projectile-detection experiment, you will get systematic errors in chi-squared
that can be huge in some cases.

3. If you are able to accurately change these experiments to the equivalent
projectile scattering range, then you can normalize them together. If you do,
read the newest Gosia manual (OP,YIEL section) to see how to calculate the
normalization constants. Even if you measure them, you will need extra factors
that are not intuitive (but easy to calculate).

Ok I will define projectile angles, that's not a big deal, although it's not so convenient. But I was really not aware that I have to calculate those normalization constants by hand. I was of the opinion as long as I have all LN=1 the experiments will automatically be cross-normalized. I found the formula on p.177 with the Rutherford cross-section and the efficiencies. My yields are corrected, so e_Ge=1, but I have to calculate those guys. I will try to use your Rachel GUI as I have installed it that might be nice.

4. Be careful if you change to the equivalent projectile scattering range. I
haven't looked at your kinematics, but this can present accuracy problems,
which are described in the scattering theory section of the new manual. If you
don't want to do this, then normalize each experiment individually by setting
LN to its own number:
...
Since gosia2 normalizes beam to target excitation, you will still have enough
sensitivity to the matrix elements, if you don't have too many fit parameters.

I am not completely sure about that. If I have the experiments independently normalized to each other, I have just one measured yield for the projectile and to reproduce it the transitional and diagonal matrix element can be varied free but they are correlated. What I want to say a smaller Q and a bigger B(E2) will give an equally good minimization I guess. So I think I should definitely calculate those normalization constants and use the procedure where all LN=1. I could also try to introduce more experiments, I think the data is rich enough.

----------------------------------------------------------------------------------

Do you see major conflicts between some of the data? For instance, do you see
that a lifetime and a fitted matrix element are in large conflict?

So far not, but obviously I did something wrong...

You could try a test fit with more parameters locked. Do you get a sensible
correlated error in the one-dimensional method, and roughly the same error in
the 2D method? I get a lot of insight by calculating yields with OP,INTI and
then putting them back in the *.yld files to do test fits. If you can't do a
fit with simulated data, then you almost certainly will have a problem with the
real data, using the same fit parameters, etc.

Can you take a few more points around the minimum in the 1-D error calculation?
I see only 1 point within the 1-sigma limits on the pdf plot, if I understand
correctly. The shape may have a smaller slope around the minimum, and going up
to chi-squared = 10 probably doesn't help you to understand the shape around
the minimum.

How are the fits terminating? Sometimes OP,MINI parameters need to be
adjusted, such as if the fit terminates because it hit chi-squared=1. and the
best fit value is a little less than 1.

I will check those things again after fixing the other stuff...

-----------------------------------------------------------------------------------

You have a weight of 10 for the spectroscopic data (lifetimes, etc.). Do you
mean to weight them so heavily? This adjustable weight can be somewhat useful,
but I don't think that there is a real mathematical justification to pick any
number other than 1. The error bars should take care of the weighting, in my
opinion; changing the weight seems no different from changing the error bars,
and that should be done individually for each measurement, if at all. If there
is no justification I haven't thought of, then increasing their weight
artificially reduces the correlated error measurement.

As I want to normalize to the target I wanted to "fix" it a little more. When I write this it sounds stupid to me, weighting is back to 1.

You could try using OP,MAP every time you use OP,MINI. This will update for
the present set of matrix elements. I have found that this can be done in the
same step, preceding OP,MINI with OP,MAP:

OP,REST
0,0
OP,MAP
OP,MINI
2100....

Uff, get an error here!

Code: Select all

ERROR- MAXIMUM SCATTERING ANGLE IS   19.96 DEGREES FOR EXPERIMENT  1
               ********* END OF EXECUTION **********


I have no idea why including the OP,MAP should give me that error, which doesn't appear at any step before. Doing the steps INTI, MAP, MINI with NTAP=3,0,4 one after the other doesn't give me that error.

It would be worth trying to run OP,CORR (with OP,INTI and NTAP = 3) at each
minimization step. Sometimes this is not too important, but it will be if the
correction factors change significantly as a function of the matrix elements.

That's what I do, in each step I first do the INTI, then the MINI command.

[hayes]

The experiments are all from the same beam run with the same energy. I don't understand completely what to put as mean energy and mean angle. Of course, they have to be in the corresponding range. The reason why I used different energies is that once somebody told me to vary them a little bit to get correction coefficients close to 1 after the integration.

If I understand you, varying the energies for that reason doesn't sound like a good idea. I'm also not sure why the correction coefficients being near 1 would be desirable. (They usually are near 1 in Gosia 1, but an overall factor should be irrelevant.) I think it would be best to use the same mean energy for each angle range in a beam run. Let me know if I'm missing something here.

I just use the ordinary mean, not weighted in any way: theta_mean = (theta_min - theta_max)/2. I used to worry about weighting somehow, according the expected Coulex probabilities, but I have tested the agreement between the fit with point calculations and the full calculations and found no conflicts.

But I was really not aware that I have to calculate those normalization constants by hand. I was of the opinion as long as I have all LN=1 the experiments will automatically be cross-normalized. I found the formula on p.177 with the Rutherford cross-section and the efficiencies. My yields are corrected, so e_Ge=1, but I have to calculate those guys. I will try to use your Rachel GUI as I have installed it that might be nice.

I really can't see why you should have to calculate these point cross sections either, since Gosia can easily do it. Maybe it was more intuitive this way to Tomasz when he wrote it. They do have to correspond to the mean scattering angle you are using in EXPT, so if you change the mean angle, update the normalization constant.

I just put Rachel 1.1.5 online, in case you have an old one. Please send feedback!

4. Be careful if you change to the equivalent projectile scattering range. I
haven't looked at your kinematics, but this can present accuracy problems,
which are described in the scattering theory section of the new manual. If you
don't want to do this, then normalize each experiment individually by setting
LN to its own number:
...
Since gosia2 normalizes beam to target excitation, you will still have enough
sensitivity to the matrix elements, if you don't have too many fit parameters.

I am not completely sure about that. If I have the experiments independently normalized to each other, I have just one measured yield for the projectile and to reproduce it the transitional and diagonal matrix element can be varied free but they are correlated. What I want to say a smaller Q and a bigger B(E2) will give an equally good minimization I guess. So I think I should definitely calculate those normalization constants and use the procedure where all LN=1. I could also try to introduce more experiments, I think the data is rich enough.

That definitely sounds right to me. It's clear to me now that you say it; maybe I didn't count your parameters correctly. Requiring the normalization will always help, anyway.

You could try using OP,MAP every time you use OP,MINI. This will update for
the present set of matrix elements. I have found that this can be done in the
same step, preceding OP,MINI with OP,MAP:

OP,REST
0,0
OP,MAP
OP,MINI
2100....

Uff, get an error here!

Code: Select all

ERROR- MAXIMUM SCATTERING ANGLE IS   19.96 DEGREES FOR EXPERIMENT  1
               ********* END OF EXECUTION **********


I have no idea why including the OP,MAP should give me that error, which doesn't appear at any step before. Doing the steps INTI, MAP, MINI with NTAP=3,0,4 one after the other doesn't give me that error.

That is very strange, but thank you for telling me. We know the source of this error (when it's not user-error), but I'm still trying to find a simple rule that will always avoid it. As Nigel explained to me, it's due to rounding of the angle transformed from one frame to another near where the Jacobian blows up, then transformed back, where it defines the other angle to be imaginary.

If it happens again, the usual way to correct it is to use the mean angle as I wrote above in each EXPT entry, and if that doesn't work, move it up or down a little. (You have to change the YNRM constants again after you find the right angle.)

Please tell me if you don't have luck using the mean scattering angle, so I can be sure to avoid it in the automated setup. So far, Rachel does well, but maybe I've just missed problem cases.

[LPGaffney]

Hi all,

I think I can confirm the behaviour seen by Christopher with regards to the mean scattering angle error in rachel/gosia. This error only occurs when OP,MAP is run at the same time as OP,MINI. If they are ran separately, i.e OP,MAP + OP,EXIT and then OP,MINI + OP,EXIT, then there is no error reported by gosia.

I don't know why this is, and whether it is trapping a real error, but I can't see it being a simple rounding error if this behaviour is being reproduced.

Cheers,
Liam