go-hep / hep

hep is the mono repository holding all of go-hep.org/x/hep packages and tools
https://go-hep.org
BSD 3-Clause "New" or "Revised" License
230 stars 35 forks source link

fastjet: rapidity is infinity for certain jets #97

Closed Bastiantheone closed 7 years ago

Bastiantheone commented 7 years ago

I noticed that for certain values the rapidity is Infinity, which cannot be plugged into the triangulation. For instance when Z and E are the same: X: -0, Y: -0, Z: 2.0866e-06, E: 2.0866e-06 Rapidity: +Inf or X: -0, Y: -0, Z: -9.05e-08, E: 9.05e-08 Rapidity: -Inf

Bastiantheone commented 7 years ago

The C++ version sets the rapidity in the following way

if (this->E() == abs(this->pz()) && _kt2 == 0) {
        // Point has infinite rapidity -- convert that into a very large
        // number, but in such a way that different 0-pt momenta will have
        // different rapidities (so as to lift the degeneracy between
        // them) [this can be relevant at parton-level]
        double MaxRapHere = MaxRap + abs(this->pz());
        if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
     } else {
       // get the rapidity in a way that's modestly insensitive to roundoff
       // error when things pz,E are large (actually the best we can do without
       // explicit knowledge of mass)
       double effective_m2 = max(0.0,m2()); // force non tachyonic mass
       double E_plus_pz    = _E + abs(_pz); // the safer of p+, p-
       // p+/p- = (p+ p-) / (p-)^2 = (kt^2+m^2)/(p-)^2
       _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
       if (_pz > 0) {_rap = - _rap;}
     }
sbinet commented 7 years ago

sounds good. I actually started to do it (then probably figured I shouldn't do it before I actually hit the issue): https://github.com/go-hep/hep/blob/master/fastjet/jet.go#L62

Bastiantheone commented 7 years ago

Oh okay I see.

It looks like the C++ version turns rap negative for different conditions in the if (pz<0) and the else(pz>0), but in the go version it is the same condition(pz>0). Is that on purpose?

sbinet commented 7 years ago

I don't know what I was thinking when I wrote that, I suspect it was initially written as you did in the PR and then factored it out of the 2 branches of the if/else, papering over the < / > difference...