Phasespace generation from decay chains

[jonas-eschle]

In order to make the package compatible with phasespace requires each Particle to know the distribution that parametrizes its mass if it is not fixed (e.g. BreitWigner, Gauss).

Where should this best be added?

[eduardo-rodrigues]

Would this not rather be a thing for a decay chain rather than a particle? Could you give an example of what and how the info would be used down the line in phasespace? Just so that I see things from the user perspective, which is crucial I reckon. Let's talk in person … :-)

[eduardo-rodrigues]

OK, down to one example. Imagine we look at B0 -> K*0 gamma:

>>> dm1 = DecayMode(0.000043300 , 'K*0 gamma', model='HELAMP', model_params=[1.0, 0.0, 1.0, 0.0])
>>> dm2 = DecayMode(0.6657, 'K+ pi-', model='VSS')
>>> dc = DecayChain('B0', {'B0':dm1, 'K*0':dm2})
>>> dc.to_dict()
{'B0': [{'bf': 4.33e-05,
   'fs': [{'K*0': [{'bf': 0.6657,
       'fs': ['K+', 'pi-'],
       'model': 'VSS',
       'model_params': ''}]},
    'gamma'],
   'model': 'HELAMP',
   'model_params': [1.0, 0.0, 1.0, 0.0]}]}

A simple way to provide you info would be via a new dedicated zfit metadata field, say

>>> dm2 = DecayMode(0.6657, 'K+ pi-', model='VSS', zfit={'model':'rel-BW', 'params':[]})
>>> dc = DecayChain('B0', {'B0':dm1, 'K*0':dm2})
>>> dc.to_dict()
{'B0': [{'bf': 4.33e-05,
   'fs': [{'K*0': [{'bf': 0.6657,
       'fs': ['K+', 'pi-'],
       'model': 'VSS',
       'model_params': '',
       'zfit': {'model': 'rel-BW', 'params': []}}]},
    'gamma'],
   'model': 'HELAMP',
   'model_params': [1.0, 0.0, 1.0, 0.0]}]}

Shall we go in this direction at least as a first attempt?

[eduardo-rodrigues]

Let me know, @mayou36, when you want to further discuss this.

[jonas-eschle]

Hi @eduardo-rodrigues, sorry for the long silence!

I tried around with some implementation. Just to be sure, the field after the particle name (the list, e.g. after 'B0') is for different possible decays? Or why is that exactly a list?

One question, technical: how to best parse it? in https://github.com/scikit-hep/decaylanguage/blob/master/decaylanguage/decay/decay.py#L564 it's a for loop and searching whether the decay happens. So there is no direct way to traverse the decay in a tree-like structure using the wrappers (DecayMode/DecayChain)?
We can traverse it with the dict, but then it's hard to get the particle name out, as there are also other keywords (e.g. model, model_params).
E.g. in dc.to_dict() how to know that the 'B0' is the particle? (in this case, it's also the mother, but more general, e.g. the 'K*0'.)

It's all surely doable somehow, but I wonder about the most elegant (and "future-proof") way

[jonas-eschle]

So I've got a student who would be interested actually to work on this (maybe, we're still figuring out), this would still be a desired feature I assume? Or did plans somehow change?

[eduardo-rodrigues]

So I've got a student who would be interested actually to work on this (maybe, we're still figuring out), this would still be a desired feature I assume? Or did plans somehow change?

Absolutely, still in planning. Thanks a lot for the ping.

That's great to hear. As you know I can't just now but let me then get back to you next week on this "thread". (In fact the comments you made above had forced me to look into some machinery a while back and I then got side tracked by a feature/bug-ish I had found. I then forgot about the starting point :D. In short, talk to you really soon.)

[simonthor]

Hello! I am the student @mayou36 mentioned in a previous message.

Currently, I have implemented a small prototype function that converts a dict in the decaylanguage format into a GenParticle instance. The function can be found in this file.
The function can be used like

gen_particle = recursively_replace(decaychain.to_dict())

I can also implement a function that directly converts a DecayChain to GenParticle, if that is preferred.

I was wondering, is the plan to convert dicts, which can contain many ways for a particle to decay, or DecayChains, which only contain one possible decay mode for each particle? Or have I misunderstood how DecayChains work? I assume that we want to convert dicts, since they represent more general decays?

[eduardo-rodrigues]

Hello @simonthor. Excellent to see work on this! I'm away on holidays. Let me get back to you upon my return ...

[simonthor]

Ok, enjoy your holidays!

I will post some more questions/problems here in the meantime:

Currently, I am using the particle package to get the mass and width of the mass distribution for a particle, i.e.,

p = Particle.from_string(name)
p.mass    # get the mass
p.width   # get the width

If there is a non-zero width, I create a mass distribution function (currently a normal distribution but this can be changed later). However, the top particle of a decay chain in phasespace must have a fixed mass, and cannot in other words have a width. This is a problem, since all particles except 6 in particle have a non-zero width.

Do you have any suggestions for how this should be solved?

Here are some of my potential solutions:

• A tolerance is set as a minimum mass width (e.g. 1e-10). All particles with a smaller mass width will be deemed as having a determined mass.
• Always use particle.mass as the determined mass and ignore the width for the top particle.
• Make it possible for the user to specify which method of the two above should be used.

I also opened an issue about this in my repo with the same information.

[simonthor]

After some discussions with @mayou36, the two questions mentioned earlier have been cleared up.
For reference, here are the conclusions:

Regarding the first question, I will implement a new object containing all possible decays as GenParticle instances. This object will be created from a dict with the decaylanguage format.

For the second question, I will implement a tolerance for the mass width that can be changed by the user. Any width smaller than the tolerance will be seen as constant.

I might post some more questions here later, once I encounter some other problems!

[eduardo-rodrigues]

Currently, I am using the particle package to get the mass and width of the mass distribution for a particle, i.e.,

p = Particle.from_string(name)
p.mass    # get the mass
p.width   # get the width

That method is legacy and the much better way is via Particle.find() for a single match or Particle.findall() if you want to find all matches.

[eduardo-rodrigues]

I might post some more questions here later, once I encounter some other problems!

Sure, feel totally free! Apologies that I'm slow as I will have the 2nd part of my holidays but after that I'm sure we can converge.

[simonthor]

That method is legacy and the much better way is via Particle.find() for a single match or Particle.findall() if you want to find all matches.

Ok, I will use find or findall then!

[simonthor]

Currently, when a particle has a mass width, I assume that the mass distribution is a normal distribution. Is there a way find out what the mass distribution is, i.e., if it is a normal distribution, BreitWigner, etc.?
Is there for example an attribute in the particle package which gives this information?
Or can I use the model parameter in decaylanguage, e.g., HELAMP, to somehow determine the mass distribution?

[eduardo-rodrigues]

Currently, when a particle has a mass width, I assume that the mass distribution is a normal distribution. Is there a way find out what the mass distribution is, i.e., if it is a normal distribution, BreitWigner, etc.?

There is no "unique solution" though amplitude analysis experts will tell you that for some popular particles it is customary to use shape "X" or Y". Most are of course motivated physically. As a default a Breit-Wigner is not a bad choice.

Is there for example an attribute in the particle package which gives this information?

No since that's not a particle property per say. Such info would go beyond what Particle is about, and more the scope of a decay related package such as this one. In fact, be aware that a decay model such as the ones in .dec files is a different thing from a modelling invariant mass distribution.

Or can I use the model parameter in decaylanguage, e.g., HELAMP, to somehow determine the mass distribution?

I don't really see how since different things. Just take the example

Decay K*+
0.6657      K0  pi+                        VSS;
0.3323      K+  pi0                        VSS;
0.0020      K+  gamma                      VSP_PWAVE;
Enddecay

where the models are different per decay mode (spins of particles are not the same).

I think the mass shape (model) is something the user will have to specify as guesses are more dangerous than anything here in the end. That's my 2 cents. But check with amplitude analysis experts if there is some way to build a "map" of best guesses that you can make available, if you really want to ...

[eduardo-rodrigues]

Hi @simonthor, @mayou36, I realise that this thread is already long and I may have missed to clarify certain points. Do let me know what you need from me at this stage, and how the API is getting along for this link between our 2 packages.

On

If there is a non-zero width, I create a mass distribution function (currently a normal distribution but this can be changed later). However, the top particle of a decay chain in phasespace must have a fixed mass, and cannot in other words have a width. This is a problem, since all particles except 6 in particle have a non-zero width.

Do you have any suggestions for how this should be solved?

I don't quite understand the issue here. can you explain? Maybe sorted in the meantime.

[eduardo-rodrigues]

Currently, I have implemented a small prototype function that converts a dict in the decaylanguage format into a GenParticle instance.
...
I can also implement a function that directly converts a DecayChain to GenParticle, if that is preferred.

You should only provide the minimal (though comprehensive for the specific goals) API. DecayChain can digest the specially-defined dicts and can also output them, so a conversion dict-GenParticle seems the right thing to me since you can then do anything. You both probably arrived at the same conclusion but let me know if there's more to it.

I was wondering, is the plan to convert dicts, which can contain many ways for a particle to decay, or DecayChains, which only contain one possible decay mode for each particle? Or have I misunderstood how DecayChains work? I assume that we want to convert dicts, since they represent more general decays?

Correct. The dict representations are the most general ones since they can describe a decay chain and a set of decay chains as typically parsed from a .dec file.

[eduardo-rodrigues]

Old question from @mayou36:

I tried around with some implementation. Just to be sure, the field after the particle name (the list, e.g. after 'B0') is for different possible decays? Or why is that exactly a list?

That's right. I tried various ways to represent the most general set of decay chains one may want and the present format seems to bo the job :-).

One question, technical: how to best parse it? in https://github.com/scikit-hep/decaylanguage/blob/master/decaylanguage/decay/decay.py#L564 it's a for loop and searching whether the decay happens. So there is no direct way to traverse the decay in a tree-like structure using the wrappers (DecayMode/DecayChain)?

There is, actually, and in fact https://github.com/scikit-hep/decaylanguage/blob/master/src/decaylanguage/decay/decay.py#L564 shows how to do it. Unless I'm misunderstanding you ... BTW, if some of the internal functions need to be more easily exposed, such as recursively_replace, I can certainly expose it via _recursively_replace; just le me know (similarly to many internal helpers available in https://github.com/scikit-hep/decaylanguage/blob/master/src/decaylanguage/dec/dec.py).

We can traverse it with the dict, but then it's hard to get the particle name out, as there are also other keywords (e.g. model, model_params).
E.g. in dc.to_dict() how to know that the 'B0' is the particle? (in this case, it's also the mother, but more general, e.g. the 'K*0'.)

This I do for the viewer, see https://github.com/scikit-hep/decaylanguage/blob/master/src/decaylanguage/decay/viewer.py. It certainly looks a but cryptic at first but it does the navigation correctly :-).

This makes me realise that for searching functionality (not that much available so far) I might reuse quite a bit of these helpers I'm referring you to!

Let's iterate as needed. Is there a little notebook you could share to see how your work is shaping up?

[simonthor]

There is no "unique solution" though amplitude analysis experts will tell you that for some popular particles it is customary to use shape "X" or Y". Most are of course motivated physically. As a default a Breit-Wigner is not a bad choice.

Ok, I might replace the normal distribution with a Breit-Wigner as a default.

I think the mass shape (model) is something the user will have to specify as guesses are more dangerous than anything here in the end.

This was something that @mayou36 also agreed with, so I will use the suggested zfit parameter that can be added to a decaylanguage dict.

I will post another comment in this issue once I have implemented the support for the zfit parameter, or if I encounter another problem. After that, this issue should mostly be resolved.

[simonthor]

If there is a non-zero width, I create a mass distribution function (currently a normal distribution but this can be changed later). However, the top particle of a decay chain in phasespace must have a fixed mass, and cannot in other words have a width. This is a problem, since all particles except 6 in particle have a non-zero width.
Do you have any suggestions for how this should be solved?

I don't quite understand the issue here. can you explain? Maybe sorted in the meantime.

This has been solved. The top particle in a phasespace GenParticle must have a constant mass, so what I do right now is to always use particle.mass and ignore the width, i.e. assume that the top particle always has a fixed mass.

[eduardo-rodrigues]

I will post another comment in this issue once I have implemented the support for the zfit parameter, or if I encounter another problem. After that, this issue should mostly be resolved.

Sounds good. Looking forward to it. (Anything like the zfit parameter I suggested at the top of this discussion?)

[simonthor]

Sounds good. Looking forward to it. (Anything like the zfit parameter I suggested at the top of this discussion?)

Yes, almost the exact way you suggested it to be implemented.

[simonthor]

Ok, the previous comment turned out to be wrong, as the implementation barely resembles the zfit parameter mentioned above. The working implementation I have added can be used as follows:

dec_dict = {'D+': [{'bf': 1,
                        'fs': ['K-', 'pi+', 'pi+',
                               {'pi0': [{'bf': 1, 'fs': ['gamma', 'gamma']}]},
                               ],
                        'model': 'PHSP', 'model_params': ''}]
                }
FullDecay.from_dict(dec_dict, {'K-': 'rel-BW', 'pi+': 'BW', 'pi0': 'gauss'})

In other words, the user passes a second dict which specifies which particle should have what type of mass function.
Is this a sufficient implementation?

[eduardo-rodrigues]

Ok, the previous comment turned out to be wrong, as the implementation barely resembles the zfit parameter mentioned above. The working implementation I have added can be used as follows:

dec_dict = {'D+': [{'bf': 1,
                        'fs': ['K-', 'pi+', 'pi+',
                               {'pi0': [{'bf': 1, 'fs': ['gamma', 'gamma']}]},
                               ],
                        'model': 'PHSP', 'model_params': ''}]
                }
FullDecay.from_dict(dec_dict, {'K-': 'rel-BW', 'pi+': 'BW', 'pi0': 'gauss'})

In other words, the user passes a second dict which specifies which particle should have what type of mass function.
Is this a sufficient implementation?

Hi @simonthor and @mayou36, thanks for the example of where you're heading.

That will work. You have, though, taken a different approach, as above I incorporate zfit specific metadata at decay mode level inside my decay decay whereas you take the approach of "composition". You loose generality - same shape for a given particle - but same space; the latter is most probably not relevant since you are not going to deal with very deeply nested chains. I confess that I'm more convinced by my suggestion, biased of course. Can you dwell on your motivations?

[simonthor]

@eduardo-rodrigues your suggested implementation seems better, since it allows for more customization, as the same particle can have different mass functions for different decay modes. To be honest, the main reason I implemented it this way was because I did not fully understand your example with the zfit parameter. More specifically, which particle does the zfit parameter correspond to? I now realize that it probably refers to the particle that decays into that decay mode, i.e., K*0 in the example at the beginning of the thread? If I have understood it correctly, your suggestion will probably be easier to implement than my current implementation, so I should be able to switch it quite quickly. Though I am also currently wrapping up my Google Summer of Code project, so unfortunately I cannot guarantee any time frame.

[eduardo-rodrigues]

Morning. Cool. But let's hear anyway from @jonas-eschle to make sure that we converge on the "best" solution to the problem for the benefit of all :-).

To be honest, the main reason I implemented it this way was because I did not fully understand your example with the zfit parameter. More specifically, which particle does the zfit parameter correspond to? I now realize that it probably refers to the particle that decays into that decay mode, i.e., K*0 in the example at the beginning of the thread?

Correct - the zfit parameter is added to DecayMode code

>>> dm2 = DecayMode(0.6657, 'K+ pi-', model='VSS', zfit={'model':'rel-BW', 'params':[]})
>>> dc = DecayChain('B0', {'B0':dm1, 'K*0':dm2})

and then becomes a "property metadata" of the particle that decays via that decay mode. This make sense. It also indeed mean that, for whatever reason, you wanted to give another K0 another mass shape, that would be possible because either the K0 would decay to a different decay mode or you would have flagged as "MyOtherK*0.

If I have understood it correctly, your suggestion will probably be easier to implement than my current implementation, so I should be able to switch it quite quickly. Though I am also currently wrapping up my Google Summer of Code project, so unfortunately I cannot guarantee any time frame.

Good luck there!

[simonthor]

Morning. Cool. But let's hear anyway from @jonas-eschle to make sure that we converge on the "best" solution to the problem for the benefit of all :-).

Ok, I will leave my implementation as it is for now then. I will discuss with Jonas on Friday, so we will probably come to a conclusion then.

I noticed one issue with your suggested implementation: how should particles withut any decays have a specified mass function?
For example, if we have the decay

{'D+': [{'bf': 1,
         'fs': ['K-', 'pi+', 'pi+',
                {'pi0': [{'bf': 1, 'fs': ['gamma', 'gamma']}]},
               ],
         'model': 'PHSP', 'model_params': ''}]
}

and want K- to have a gaussian mass distribution. How should this be specified? Should it e.g., be specified like this:

{'D+': [{'bf': 1,
         'fs': [
                {'K-': [{'zfit': 'gauss'}]},
                'pi+', 'pi+',
                {'pi0': [{'bf': 1, 'fs': ['gamma', 'gamma']}]},
               ],
         'model': 'PHSP', 'model_params': ''}]
}

This would introduce particles which are dicts that do not have fs and bf fields. Would this break backwards-compatibility for decaylanguage?

[eduardo-rodrigues]

I noticed one issue with your suggested implementation: how should particles withut any decays have a specified mass function?

My point is that a particle only needs a mass function if it decays. A proton never has/needs a mass function, for example :-).

This would introduce particles which are dicts that do not have fs and bf fields. Would this break backwards-compatibility for decaylanguage?

Yes in a way, but you can in fact use the from_dict constructor to check that straightforwardly with some design proposals, see https://github.com/scikit-hep/decaylanguage/blob/master/src/decaylanguage/decay/decay.py#L414,

As you can see the fs field expects a list of strings, say 'fs': ['K+', 'pi-'], and an element in the list that is a dict is interpreted as a decay mode. I actually did not try an empty decay mode recently, meaning one with no given final-state particles, though I'm pretty sure that I test such possibilities or at least looked into that in the past, and the code probably digest such cases fine. Can you give it a try?

If what you suggest is really needed then we would have to think harder if the present code fails to parse. If you get stuck I can try and help next week, BTW ...

[simonthor]

I talked with Jonas on Friday (i.e., yesterday) and he agreed with you that only particles that decay need a mass function; all stable particles should have a constant mass.
I have now implemented this zfit parameter into the converter. This can be used like so:

dec_dict = {'D+': [{'bf': 1,
         'fs': ['K-', 'pi+', 'pi+',
                {'pi0': [{'bf': 1, 'fs': ['gamma', 'gamma']}]},
               ],
         'model': 'PHSP', 'model_params': '',
         'zfit': 'rel-BW'}]
}

decay = FullDecay.from_dict(dec_dict)

The last line will create a FullDecay object which is a small wrapper around a list of GenParticle instances and their probablilities. One can then call decay.generate, which takes the same arguments as GenParticle.generate and will call that method under the hood.

What I have left now is to "only" implement some mass functions, such as the relativistic Breit-Wigner (right now I only have a Gaussian implemented). I plan on using zfit (based on a suggestion from Jonas) to create these mass functions. Do you have any suggestions for which mass functions I should implement? By the way, it is already possible for a user to add additional mass functions, besides the common ones I plan to implement.

Are there any more features that you would like to see implemented? Otherwise, implementing some specific mass functions will be the last feature that needs to be added. Then this issue should mostly be solved.

[simonthor]

To find the mass of a particle, I currently use Particle.findall in a somewhat clunky way, since Particle.find was deprecated (code copied from my repo):

def _find_single_particle(name: str) -> Particle:
    particles = Particle.findall(lambda p: p.name == name)
    if len(particles) != 1:
        raise ValueError(f'The given name "{name}" of the particle was not valid.'
                         'The input dict probably had an invalid particle name in it.')
    return particles[0]

Can I simply use Particle.from_evtgen_name instead? Are EvtGen names always used for particle names in decaylanguage dicts? If not, is there a different method that you suggest me to use, tbat has a cleaner syntax than what I do now?
If this is better suited as an issue on the particle package, I can of course move this comment there.

[jonas-eschle]

I have now implemented this zfit parameter into the converter. This can be used like so:

That seems very reasonable to me so far. I think a relativistic BW, a normal BW and a Gaussian seem to be useful.

@eduardo-rodrigues I was wondering a bit on the logistics of this, as it uses multiple packages and therefore does not quite fit into one single package (we could add it directly to phasespace, that is maybe the best bet...) and seems to be nearly a standalone package. Or where would you see this fit?

Furthermore, it made me think about a RapidSim version. We have the phasespace, the dec parser, so we need some efficiencies, smearing and we would be able to do rapid sim. Do you think this is another possible project (for Simon or anyone else)? I think it could be quite useful

[eduardo-rodrigues]

I talked with Jonas on Friday (i.e., yesterday) and he agreed with you that only particles that decay need a mass function; all stable particles should have a constant mass.

Good. Reassuring ;-).

I have now implemented this zfit parameter into the converter. This can be used like so:

dec_dict = {'D+': [{'bf': 1,
         'fs': ['K-', 'pi+', 'pi+',
                {'pi0': [{'bf': 1, 'fs': ['gamma', 'gamma']}]},
               ],
         'model': 'PHSP', 'model_params': '',
         'zfit': 'rel-BW'}]
}

If I understand you correctly the above requires no coding on your side at all. Indeed

In [1]: from decaylanguage import *

In [2]: dm1 = DecayMode(1, 'K- pi+ pi- pi0', model='PHSP', zfit='rel-BW')

In [3]: dm2 = DecayMode(1, 'gamma gamma')

In [4]: dc = DecayChain('D+', {'D+':dm1, 'pi0':dm2})

In [5]: dc.to_dict()
Out[5]:
{'D+': [{'bf': 1,
   'fs': ['K-',
    'pi+',
    'pi-',
    {'pi0': [{'bf': 1,
       'fs': ['gamma', 'gamma'],
       'model': '',
       'model_params': ''}]}],
   'model': 'PHSP',
   'model_params': '',
   'zfit': 'rel-BW'}]}

which is exactly what you wrote. The benefit is that my construction is trivial to read by anyone, even a reasonably newbie.

decay = FullDecay.from_dict(dec_dict)

The last line will create a `FullDecay` object which is a small wrapper around a list of GenParticle instances and their probablilities. One can then call `decay.generate`, which takes the same arguments as `GenParticle.generate` and will call that method under the hood.

Here I cannot judge too well but if this is the way you want to have the link/bridge between the 2 packages, then fine.

What I have left now is to "only" implement some mass functions, such as the relativistic Breit-Wigner (right now I only have a Gaussian implemented). I plan on using zfit (based on a suggestion from Jonas) to create these mass functions. Do you have any suggestions for which mass functions I should implement? By the way, it is already possible for a user to add additional mass functions, besides the common ones I plan to implement.

Not apart from the ones you mention. I guess your intention is to provide a way for the user to pass a user-defined function where you would map a name MyCoolMassFunction to the actual function? Thinking out loud.

Are there any more features that you would like to see implemented? Otherwise, implementing some specific mass functions will be the last feature that needs to be added. Then this issue should mostly be solved.

From the above I wonder if you would like FullDecay.from_DecayChain(dc). But I guess not since FullDecay.from_dict(my_DecayChain_instance.to_dict()) will do trivially. This is in fact best since minimising the API - I also only make the viewer understand dicts for the same reason.

[eduardo-rodrigues]

To find the mass of a particle, I currently use Particle.findall in a somewhat clunky way, since Particle.find was deprecated (code copied from my repo):

def _find_single_particle(name: str) -> Particle:
    particles = Particle.findall(lambda p: p.name == name)
    if len(particles) != 1:
        raise ValueError(f'The given name "{name}" of the particle was not valid.'
                         'The input dict probably had an invalid particle name in it.')
    return particles[0]

Can I simply use Particle.from_evtgen_name instead?

For sure. This is very handy when dealing with decay files since they use EvtGen names rather than the usual PDG ones (I know, these matters are annoying, but not much one can do at this point in history). The DecayChain class and related do not care except for the charge conjugation method, see the docs there ...
To be comprehensive: there might be subtleties - I deal with a bunch in this package if you go through the code - but better not complicate it too much. We can talk about those if needed after you're basically done. For example, if a user defines a decay with a PDG name, or defines "MyK+" and you need to know this is a "K+" ...

Are EvtGen names always used for particle names in decaylanguage dicts? If not, is there a different method that you suggest me to use, tbat has a cleaner syntax than what I do now?

So by default DecayLanguage dicts will deal with EvtGen names. For sure if they are constructed from .dec files. This being said, the dicts can equally deal with PDG names and only some methods need to know about that when called, see e.g. https://github.com/scikit-hep/decaylanguage/blob/master/src/decaylanguage/decay/decay.py#L76. (Dunno how much you want my replies to be verbose.)

If this is better suited as an issue on the particle package, I can of course move this comment there.

Happy to discuss anything else there if better. Till now it seems good to have all the discussion of this bridging between the 2 packages in this one place.

[eduardo-rodrigues]

@eduardo-rodrigues I was wondering a bit on the logistics of this, as it uses multiple packages and therefore does not quite fit into one single package (we could add it directly to phasespace, that is maybe the best bet...) and seems to be nearly a standalone package. Or where would you see this fit?

Having it in phasespace also seems the natural place for me. In a way what we are doing here is not too different from the bridges we have designed between ROOT, uproot, Hist and boost-histogram, etc., which acknowledge the relevant representations in the "PyHEP ecosystem". With this we are effectively accepting that this ecosystem deals with particles and decays with Particle and DecayLanguage, and phasespace can generate based on a decay construction provided by DecayLanguage. As for dependencies it is also not different from, say, some plotting functionality in Hist requiring matplotlib only if you want to make plots as well, see https://github.com/scikit-hep/hist/blob/main/src/hist/plot.py and the verification at https://github.com/scikit-hep/hist/blob/main/setup.cfg#L43-L47 that no plotting library is required by default. In this case you would only require DecayLanguage for the functionality discussed here. Agreed?

Furthermore, it made me think about a RapidSim version. We have the phasespace, the dec parser, so we need some efficiencies, smearing and we would be able to do rapid sim. Do you think this is another possible project (for Simon or anyone else)? I think it could be quite useful

This could be a good project but one would need to think carefully on how to make the thing at least as good as RapidSim: functionality but also installation easiness and package dependencies. Maybe something for a dedicated discussion on the PyHEP channel?

[eduardo-rodrigues]

I hope I have by now replied to all outstanding points. Let me know. Seems we're not far from seeing some code in action 👍!

[eduardo-rodrigues]

Further on design. The present example in the README (BTW, I corrected a trivial typo, see zfit/phasespace#61) goes roughly as follows:

PION_MASS = 139.57018    
# ...
pion = GenParticle('pi-', PION_MASS)
kaon = GenParticle('K+', KAON_MASS)
kstar = GenParticle('K*', KSTARZ_MASS).set_children(kaon, pion)
gamma = GenParticle('gamma', 0)
bz = GenParticle('B0', B0_MASS).set_children(kstar, gamma)
weights, particles = bz.generate(n_events=1000)

For the functionality exploiting DecayLanguage a neat API, IMO, would be

dm_kstarz = DecayMode(1, 'K+ pi-', zfit='rel-BW')
dm_bz = DecayMode(1, 'K*0 gamma')
dc = DecayChain('B0', {'B0': dm_bz, 'K*0': dm_kstarz})
bz = GenParticle.from_decaylanguage(dc)
# bz = GenParticle.from_decaylanguage(dc.to_dict()) is also possible but seems more convoluted
# since the navigation via DecayChain seems easier and relating well to GenParticle.set_children() and related. Right?
weights, particles = bz.generate(n_events=1000)

This suggested class method has the advantage that your package dependency on DecayLanguage is only via here, so not for the whole package.

Following up from here you could then also have a method GenParticle.to_decaylanguage() that would allow you to visualise what is defined with DecayChainViewer(bz.to_decaylanguage().to_dict()).

Still thinking out loud: I have not checked whether one can easily get the used particle masses from GenParticle but that would be a useful and meaningful thing to have when using DecayLanguage. Indeed, unlike what happens with the README example above, where masses are explicitly passed in, for the "bridge functionality" the masses are taken from the Particle package, hence the user does not see upon creation of the GenParticle instance what masses are going to be used. Maybe worth making an issue in phasespace?

[simonthor]

Not apart from the ones you mention. I guess your intention is to provide a way for the user to pass a user-defined function where you would map a name MyCoolMassFunction to the actual function? Thinking out loud.

Yes, exactly. One can pass an additional dict to the constructor like so:

FullDecay.from_dict(dec_dict, {'cool func': MyCoolFunction})

which will use MyCoolFunction as the mass function for all decay modes with zfit: 'cool func'.

From the above I wonder if you would like FullDecay.from_DecayChain(dc). But I guess not since FullDecay.from_dict(my_DecayChain_instance.to_dict()) will do trivially. This is in fact best since minimising the API - I also only make the viewer understand dicts for the same reason.

If from_DecayChain should be implemented, it is probably better to add that functionality as a class method to GenParticle (as you mentioned in another comment). Adding that function to FullDecay is likely not necessary.

Still thinking out loud: I have not checked whether one can easily get the used particle masses from GenParticle but that would be a useful and meaningful thing to have when using DecayLanguage. Indeed, unlike what happens with the README example above, where masses are explicitly passed in, for the "bridge functionality" the masses are taken from the Particle package, hence the user does not see upon creation of the GenParticle instance what masses are going to be used. Maybe worth making an issue in phasespace?

For particles with constant mass, it is possible to get the mass with GenParticle.get_mass(). However, for particles with variable mass, this function returns a random sample from the mass distribution (if I have understood it correctly). I do not think a user can get the underlying mass function, i.e., if it is a normal distribution and what its mu and sigma is. Maybe @jonas-eschle has a better answer to this question?

[eduardo-rodrigues]

Two quick questions:

1. What exactly is the scope of your project? Within GSoC or ...? You are a student at ...? Just being curious.
2. What is your timescale to create a first PR with the design status? Looking forward to it.

[simonthor]

1. What exactly is the scope of your project? Within GSoC or ...? You are a student at ...? Just being curious.

I work on this during my free time. I am a student from Sweden (KTH Royal Institute of Technology to be specific) that aspire to become a particle physicist and I am interested in contributing to Scikit-HEP. Sorry, I should have introduced myself much earlier!

2. What is your timescale to create a first PR with the design status? Looking forward to it.

I have implemented most of what we have discussed, besides the relativistic Breit-Wigner distribution, which will likely be added to zfit-physics as a separate PR. All code I have written can be found in my repo (thought I sent the link before but it might very well have been lost in all the messages).
Hopefully, I will mostly be done with the relativistic Breit-Wigner zfit PDF by the end of this week, and I can then open a PR to phasespace.

[simonthor]

I have now added relativistic Breit-Wigner to zfit-physics and use it as the default mass function for the decaylanguage to phasespace conversion. The only task remaining now should be to merge the code into phasespace.

[jonas-eschle]

Good to hear that it works! So for that, what do you think @eduardo-rodrigues , as it is it greatly extends the capability of phasespace and is a real generator that is also capable of generating from multiple decay branches. So I see it as a very possible future to extend the functionality even further (which is difficult if it is in phasespace).

I think I would tend to make it a standalone repo, we can very well add it to the zfit org for the moment (it may not be fitted, yet, for Scikit-HEP, or what do you think?).

@simonthor what do you think? We could keep it a standalone repo, which allows to have it extendible beyond what phasespace can do.

[simonthor]

I agree that keeping it as a standgalone repo would make more sense, also since it has many dependencies that phasespace currently does not have.

[eduardo-rodrigues]

Hi folks, great news. See my comments and suggestions above on the to_decaylanguage and from_... with no need for extra dependencies, which is a "technique" that uproot, hist and many others do. At this stage, and given the small amount of code involved, I think there is more to gain from augmenting the present packages than from a separate package. My 2 cents. Happy to talk also live in Zoom if that helps, as it seems we're discussion a lot (for good reasons).

[eduardo-rodrigues]

Hi both, any follow-ups?

[simonthor]

Hello, apologies for the delayed reply.
A Zoom meeting sounds like a good idea to me. I emailed Jonas ~2 days ago but have not received a reply yet, so he might be busy right now.

[eduardo-rodrigues]

OK, fair enoug. Let's then wait.

[eduardo-rodrigues]

Closing since zfit/phasespace#63 provides a first version of the functionality.