Page MenuHomeHEPForge

Implementation of higher-spin K-matrix
ClosedPublic

Authored by johndan on Nov 26 2020, 12:36 PM.

Details

Summary

This differential contains commits intended to allow use of K-matrix implementation in non S-wave scenarios.
(it builds on D39)

Test Plan

Trial fits to a B+ 3-body decay

Diff Detail

Repository
rLAURA laura
Branch
johndan-KMatrix_D0K_higherSpin
Lint
No Lint Coverage
Unit
No Test Coverage
Build Status
Buildable 123
Build 123: arc lint + arc unit

Event Timeline

There are a very large number of changes, so older changes are hidden. Show Older Changes
  • Add K-matrix name to gCoupling parameters to avoid clashes in models with >1 K-matrix
  • Resolve merge conflict
  • Allow other parameters of the SVP part to be floated
  • Fix merge conflict

Thanks, @jback. A few answers to your responses, with just one that doesn't (yet) make sense to me...

src/LauKMatrixProdPole.cc
87

On the first point, I don't think it matters if other channels have L>0 does it? Why do we need to track it?

On the second, I think you're spot on - the resBWFactor has been dropped as you say.

src/LauKMatrixPropagator.cc
97

Hmmm not sure. If a = 3 then calcGamma gives:
pow(q,L) / pow( q*q + gamAInvRadSq_[iCh] , L_[iCh]/2. );
where gamAInvRadSq_ = a/(R*R)

i.e. the return is:

L = 2 : sqrt[ (qR)^2L / ( (qR)^2 + 3 )^L ]
L = 2 : sqrt[ (qR)^4 / ( (qR)^2 + 3 )^2 ]
L = 2 : sqrt[ (qR)^4 / ( (qR)^4 + 6 (qR)^2 + 9 ) ]

But looking at equation 98 in the Laura++ paper I would expect the coefficient of (qR)^2 in the denominator to be 3 not 6. What am I missing?

johndan marked an inline comment as done.
  • Remove duplicate par_mSq0_
src/LauKMatrixProdPole.cc
87

Yes, we don't need to track this since we have changed how the angular dependence is included.

src/LauKMatrixPropagator.cc
97

You are not missing anything; it's just that we can't exactly reproduce the z^4 + 3z^2 + 9 factor using a perfect square product between two D terms without introducing complex terms. Remember, this barrier term is an ad-hoc "best guess" at what the momentum dependence should be. The important thing is that we have a quadratic dependence on z^2, and we can change "a" and "R" at initialisation for systematic variations.

johndan added inline comments.
src/LauKMatrixPropagator.cc
97

OK! Agreed, the dependence is what matters.

johndan marked an inline comment as done.
  • Use default values at declaration in LauAbsResonance to significantly simplify the constructors
  • Protect all parameters of calls to resAmp for prodpole / svp with const
  • Clean up K-matrix propagator default values. Add lots of const protection to function calls. Noexcept throughout

The exact amplitude has been reproduced with another fitter. Anything left for me to do before we land?

  • Movement of one comment
  • Allow floating of s0scatt and mSq0
  • Cherry pick commit from 19th Jan that moved enum class KMatrixChannels and resolve associated conflicts
  • Turn off auto-verbose-printout in propagator
  • Typo in LauKMatrixPropagator.cc
  • Updates to K-matrix propagator to use refined barrier factor
  • Clean up indentation
  • Parallel updates to KMatrixPropagator.hh
  • Cleanup indentation
  • Revert "Allow slope of exponential non-resonant shape to vary positive (for low-stats fits)"
  • rebase conflict
  • rebase conflict
  • D*K system should be in P-wave\!
  • Changes after round of code-review with @jback
  • Remove whitespace
  • Remove whitespace
  • Remove extraneous comment
  • Fix compilation errors
  • rebase conflict
  • Simplify setting of denominator parameter in barrier factor
  • Cache the a/R^2 calculation
  • Neater access to D*0 mass information from LauResonanceMaker - thanks to Tom
  • Add definition of angular momentum term for spin-2 (max spin) to K-matrix
  • Remove done TODO comment and make const several function arguments
  • rebase conflict
  • Add barrier factors to numerator of K-matrix amplitude in ProdPole and ProdSVP, which need to know the spin of the K-matrix using the new getL() function
  • Fix typo in ProdPole.cc and ProdSVP.cc
  • - Changes to allow different orbital angular momenta for different channels inside the K-matrix
  • Correctly initialise a_ vector. Remove calculation of angular factor from K-matrix propagator. Remove any need to update the propagator with the LauKinematics object
  • Remove angular parts of K-matrix header, and add extra function to indicate the final-state channel of the various coupled channels
  • Update ProdPole and SVP to use resAmp (and leave amplitude to LauAbsResonance)
  • Update ProdPole and SVP headers
  • rebase conflict
  • Changes to LauResonanceMaker to allow retrieval of parent Blatt-Weisskopf factor without using resInfo that is not defined for K-matrix
  • Changes to LauBlattWeisskopfFactor to accommodate tangential use by K-matrix to get parent factor without access to a resInfo
  • Add K-matrix name to gCoupling parameters to avoid clashes in models with >1 K-matrix
  • Resolve merge conflict
  • Allow other parameters of the SVP part to be floated
  • Fix merge conflict
  • rebase conflict
  • Remove duplicate par_mSq0_
  • rebase conflict
  • Protect all parameters of calls to resAmp for prodpole / svp with const
  • resolve conflict
  • Missing ;
  • Movement of one comment and correct access to one floating parameter
  • Remove unnecessary whitespace
  • Allow floating of s0scatt and mSq0
  • Cherry pick commit from 19th Jan that moved enum class KMatrixChannels and resolve associated conflicts
  • Turn off auto-verbose-printout in propagator
  • Move a comment - trivial
  • Increase precision of verbose printout
  • Clean whitespace
  • minor change cherry-picking earlier BelleNR commit
tlatham requested changes to this revision.Feb 3 2021, 12:07 AM

Thanks for all this @johndan it looks really nice. I've not really gone through the physics since I think you and @jback have gone through that thoroughly already. So these comments are all more technical.
There are various inline comments (mostly minor). Plus, here are a few general comments:

  • Do you have an example of using this generalisation of the K-matrix that could go into the examples dir?
  • Please remove the noexcept specifications on functions
    • I would prefer to do a systematic sweep through the package as part of T34 and I still need to fully firm up my understanding of when it is necessary/helpful to do this and what are the prerequisites for it to be valid to do so
    • Having said that, I'm pretty certain that some of those that you've specified noexcept are actually potentially throwing because they call functions that are themselves potentially throwing, e.g. TMatrixD::operator[]
  • Please update the release notes
inc/LauKMatrixProdPole.hh
79–82

The Doxygen for getFloatingParameters() needs to be restored.

84–85

The first argument here is the mass not the mass squared.

inc/LauKMatrixProdSVP.hh
78–81

The Doxygen for getFloatingParameters() needs to be restored.

84–85

The first argument here is the mass not the mass squared.

inc/LauKMatrixPropagator.hh
514–515

Use uniform initialisation here for consistency.

551–585

These can all be made const I think

594–595

This can be static constexpr I think.

src/LauKMatrixProdPole.cc
150–168

Braces have gone missing

src/LauKMatrixPropagator.cc
114–116

I don't think that this loop does anything because the values of a_ should already have been set and then used to set gamAInvRadSq_ in setParameters above. If they were used anywhere else, this would be bad because they're being overwritten. But actually they don't seem to be used anywhere else, I'm not sure that these even need to be a member variable, they could be local to setParameters.
I would say that they should be all initialised to 0 at the start of that function, then there should be a loop like this one in storeOrbitalAngularMomenta but that should only be used if storeBarrierFactorParameter has not been called already - so need some sort of flag for that.

312–315

Need to mention the angular momentum and the barrier factor parameters here too.

442–454

These two vectors will have size 2*nChannels_ after this.
I would suggest this formulation instead:

L_.clear();
L_.assign( nChannels_, 0 );

and similarly for radii_.

456–458

As discussed above, probably this shouldn't be a member variable. However, if it is decided to keep it as one then can also use assign here.

552

Looks like this change from D39 got reverted by mistake here.

598

Another accidental reversion

611–644

Not sure why this wording has changed

632

Should be atoi here rather than atof

1039–1061

Might be nicer to convert this to a switch - that then has the compile-time check that there is a case for each state of the enum.

This revision now requires changes to proceed.Feb 3 2021, 12:07 AM
johndan marked 20 inline comments as done.
  • Remove noexcept
  • Update release notes
  • Restore getFloatingParameters doxygen comment
  • Fix LauKMatrixProdSVP whitespace issues
  • Fix parameter naming in resAmp function declaration for ProdPole and SVP
  • Further whitespace fix to SVP
  • Uniform default variable
  • Add const variables
  • Make verbose flag static constexpr (nice, thanks Tom\!)
  • Add braces for Tom
  • Move setting of propagator denominator parameter to setParameters(...) routine
  • Correct a initialisation and move it to a local variable in storeParameters()
  • Update setParameters(...) text to describe new non-zero-spin-related input lines in configuration file for K-matrix
  • Simplify default assignation for L_ and radii_ vectors
  • Reinstate D39 name for pole mass sq param
  • Reinstate D39 name for scatt constant param
  • Reverse change to error message
  • Use atoi for angular momenta not atof
  • Add switch on phasespaceindex inside calcGamma(...) to catch undefined calls to that function

Do you have an example of using this generalisation of the K-matrix that could go into the examples dir?

Not yet, but very happy to add one when this implementation has been tested a bit more in my fits. Can I do this later, or shall I pause D41 landing until I've gained a bit more experience with the K-matrix in D0 K+?

inc/LauKMatrixProdPole.hh
84–85

Yuk. This was correct in the implementation. Is there no compile-time check that parameter names in the declared function match those of the implementation?

src/LauKMatrixPropagator.cc
611–644

that's an error, thanks.

  • Remove noexcept
  • Merge release notes
  • Restore getFloatingParameters doxygen comment
  • Fix LauKMatrixProdSVP whitespace issues
  • Fix parameter naming in resAmp function declaration for ProdPole and SVP
  • Further whitespace fix to SVP
  • Uniform default variable
  • Add const variables
  • Make verbose flag static constexpr (nice, thanks Tom\!)
  • Add braces for Tom
  • Move setting of propagator denominator parameter to setParameters(...) routine
  • Correct a initialisation and move it to a local variable in storeParameters()
  • Update setParameters(...) text to describe new non-zero-spin-related input lines in configuration file for K-matrix
  • Simplify default assignation for L_ and radii_ vectors
  • Reinstate D39 name for pole mass sq param
  • Reinstate D39 name for scatt constant param
  • Reverse change to error message
  • Use atoi for angular momenta not atof
  • Add switch on phasespaceindex inside calcGamma(...) to catch undefined calls to that function
  • Fix 5 compile-time syntax errors
  • Return propagator pointer in LauIsobarDynamics::defineKMatrixPropagator
  • Change proposed by Tom to ensure that K-matrix is always treated as non-narrow for integration
  • Fix: set spin type correctly to Legendre for the K-matrix
  • Reset haveCalled... variable in initialiseMatrices()

All looks good to me thanks @johndan. If @jback can do a final check over as well, since there have been a few changes, then hopefully we can get this landed soon.

inc/LauKMatrixProdPole.hh
84–85

No, only the types are checked by the compiler

src/LauKMatrixPropagator.cc
973

We need to change this to find gamma for all channels, not just DK. I think we can do:

q = sqrt(s)*mag(rho)/2

where rho is the (complex) phase space factor depending on the channel, similar to what is done in calcRhoMatrix().

We should also probably move the if statements within calcRhoMatrix() to another function like getRho(s, channel), which can also be called here to avoid code repetition.

  • Factorise retrieval of the phase-space density and use to obtain barrier-factor matrix for all channels
johndan marked 3 inline comments as done.
  • Compilation fixes
src/LauKMatrixPropagator.cc
973

I think all bases were covered in that calcGamma(...) was not called (by calcGammaMatrix(...)) for any case other than DK because only the two DK channels had L!=0 and the default 'gamma' was trivially returned in other cases. But I see the general usefulness of what you propose, and that it nicely avoids redetermining the break-up momentum in calcGamma so I've implemented that as you proposed.

Thanks!

Thanks @johndan for implementing my suggestion.

Code looks ready, apart from one very minor issue:

Please remove "DK P-wave" from the comment "Calculate the DK P-wave gamma angular-momentum barrier" for the calcGamma() function, both in the source and header file of the propagator class, since the method is now general.

I'll then accept this revision.

inc/LauKMatrixPropagator.hh
269

Please remove "DK P-wave" in the comment here, since method is general.

Code can be merged, although it would be good if the "DK P-wave" comment for calcGamma in LauKMatrixPropagator.hh is removed before this revision is landed.

Thanks very much @johndan for implementing this code. A great advancement for the K-matrix amplitude method in Laura++.

This revision is now accepted and ready to land.Feb 10 2021, 9:18 PM
This revision was automatically updated to reflect the committed changes.

Just to say thanks very much to @jback and @tlatham for taking time to do a really thorough review. Much appreciated!