Particle polarization in events

[kilian]

We have to fix particle (tau) polarization conventions in event output. Also, there is more polarization info than just helicity projections.

We need handles for the polarization info that is inside Whizard, but can't be currently written out.

On the way to enabling this feature, I run into the deficiencies of Whizard's treatment of individual particle polarization by state matrices, which are fragile due to their internal pointers. I'll solve this by switching to an internal representation with allocatables. (This issue has bitten me before.) Then proceed with the I/O handles.

[Bijan Chokoufe Nejad]

Are you working on this right now? We discussed Akiya's patch yesterday and I would include his changes with the modifications that we discussed. But this will likely conflict with the changes that you are planing for the polarization. Or maybe not. He only puts the information in the helicity_t type of the particle_t but the hepeup output has to be modified as well. So should I put the changes in or wait until you are done?

[kilian]

OK, if you do have the patch, just commit it.

The changes that I'm considering would affect the polarization_t component. So that's not a problem.

I'm just slightly concerned that helicity_t is one of the fundamental types that go into state matrix, interaction, evaluator and essentially anything. Adding something inside would probably not be visible anywhere else, but one shouldn't expect the addition to be respected in raw event files unless explicitly coded. Also, if this is visible in helicity output, it will modify all MD5 sums.

[Bijan Chokoufe Nejad]

There was one issue with different conventions concerning what should be written into SPINUP. What you currently do in hepeup_from_particle_set is write in the scalar product of the three momentum of the mother with spin vector, which seems to be the LHA convention.

But this subroutine is also called in shower_pythia6_import_particle_set where the Fortran Tauola seems to expect just the helicity in SPINUP. At first we thought about a flag ps_use_tauola to enforce this second convention but it might be better to add an optional argument in the call from the shower module to demand this. Any opinion on this?

[Bijan Chokoufe Nejad]

Another change from Akiya that you might object to: In hepevt_set_particle you seem to set the spinlh 3-vector to [r, theta, phi]. He seems to expect

-          if (pol_status == PRT_DEFINITE_HELICITY) then
-             select case (hel)
-             case (1:)
-                spinlh(:,i) = [one, zero, zero]
-             case (:-1)
-                spinlh(:,i) = [one, PI, zero]
-             case (0)
-                spinlh(:,i) = [one, PI/2, zero]
-             end select
-          else
-             spinlh(:,i) = [zero, zero, zero]
-          end if
+          spinlh(:,i) = [zero, zero, zero]
+          spinlh(3,i) = hel

So just writing the helicity in the 3rd component that is interpreted as z component. In the StdHep? Manual 5.06.01, the only comment to spinlh is "spin information". Do you know any source that would vouch for a r theta phi interpretation? Otherwise I would apply the pragmatic solution.

[kilian]

This is precisely the issue that we discussed. Since I wasn't sure about the conventions, I chose one. Akiya just told us that I chose the wrong one. So correct this, unless you see a conflicting use of the conventions elsewhere. (One might also invent a parameter for switching conventions.)

In any case, please write exhaustive tests (unit or functional) for the corrected and any new functionality. Otherwise it becomes too easy to screw this up by any future change.

[Bijan Chokoufe Nejad]

OK very good. I will use polarization_get_axis (pol) for the general case then. What is not quite clear to me is when PRT_DEFINITE_HELICITY and when PRT_GENERIC_POLARIZATION is set. I found out that it's connected to FM_SELECT_HELICTY and FM_FACTOR_HELICITY, respectively. Furthermore FM_SELECT_HELICITY is set when we demand ?polarized_events. But how is the generic case activated?

[kilian]

Exactly what I want to do. Everything but PRT_SELECT_HELICITY is broken when viewed from the outside. Trying to enable the other cases, I ran into memory issues with the polarization state matrix. I'd like to fix that now.

[Bijan Chokoufe Nejad]

I see. I have Akiyas changes, which I think are necessary, now in Whizard.

But the only thing I can test right now, is the polarization in the output of a particle that is produced by Whizard. To check the whole chain, I need the tauola folder.

It would be good to have the phone meeting wether we can include this in the short term.

[kilian]

Just to note that I'm now going to refactor parts of the polarizations module, in particles.nw.

[Juergen Reuter]

Intel compiler lets commands_26 (alternate setups) fail now, probably a workaround for allocatables is necessary there:

$ diff err-output/commands_26.out ref-output/commands_26.ref 
34a35,60
>  Alternative setups:
> ------------------------------------------------------------------------
>  Commands:
> ------------------------------------------------------------------------
> +  SEQUENCE    <local_command_list>  =  <cmd_num>
> +  SEQUENCE    <cmd_num>  =  <var_name> '=' <expr>
> |  +  IDENTIFIER  <var_name>  = i
> |  +  KEYWORD     '='  = [keyword] =
> |  +  SEQUENCE    <expr>  =  <term>
> |  |  +  SEQUENCE    <term>  =  <factor>
> |  |  |  +  SEQUENCE    <factor>  =  <integer_value>
> |  |  |  |  +  SEQUENCE    <integer_value>  =  <integer_literal>
> |  |  |  |  |  +  INTEGER     <integer_literal>  = 1
> ------------------------------------------------------------------------
>  Commands:
> ------------------------------------------------------------------------
> +  SEQUENCE    <local_command_list>  =  <cmd_num>
> +  SEQUENCE    <cmd_num>  =  <var_name> '=' <expr>
> |  +  IDENTIFIER  <var_name>  = i
> |  +  KEYWORD     '='  = [keyword] =
> |  +  SEQUENCE    <expr>  =  <term>
> |  |  +  SEQUENCE    <term>  =  <factor>
> |  |  |  +  SEQUENCE    <factor>  =  <integer_value>
> |  |  |  |  +  SEQUENCE    <integer_value>  =  <integer_literal>
> |  |  |  |  |  +  INTEGER     <integer_literal>  = 2
> ========================================================================

[Juergen Reuter]

As discussed with WK this was an accidental change in whizard.nw due to a local workaround for a compiler bug. This has been reverted in r7458.

[kilian]

Implementation is under way.

[kilian]

Status update: the new internal implementation via generalized Bloch vectors is finished and consistent for all supported spins (up to 2). I'm currently checking the unit tests - some interface problems are not yet fixed, so be patient.

Overall this turned out more cumbersome to do correctly than anticipated, because there are so many special cases. The new implementation is robust (no finalizer needed, intrinsic assignment works) and even has some new physics content that we might exploit in the future.

[kilian]

New implementation is enabled in r7491.

• details in modules su_algebra and bloch_vectors
• internals hidden in module polarizations, interface essentially unchanged
• no finalizer required for polarization, assignment should work correctly
• polarization iterator and method to_state ensure compatibility with previous implementation
• some reordering of helicities in unit tests, no visible effect at higher level.

The implementation has been tested with nagfor 6 and gfortran 4.8, so far.

As soon as this works with all configurations, I can look at the factorization options.

[kilian]

To be cont'd after dealing with model scheme issues.

[Juergen Reuter]

This has been stalled for a while now. What exactly is the status hereof?

[Juergen Reuter]

Milestone renamed

[Juergen Reuter]

Milestone renamed

[Juergen Reuter]

Milestone renamed

[Juergen Reuter]

Milestone renamed

[Juergen Reuter]

Transferred to Gitlab Issue 306