[Rivet] thrust

[Andrii Verbytskyi]

Dear Andy,
my statement is that 

1)
a)
"// This function implements the iterative algorithm as described in the
// Pythia manual. We take eight (four) different starting vectors
// constructed from the four (three) leading particles to make sure that
// we don't find a local maximum.

"
is not completely correct, as in some cases Pythia algorithm, as well as
other implementations arrange objects by p_T and not by |p|. That
affects some small fraction of events, and has nothing to do with
transverse thrust.

b)
"I checked the header and documentation and 
> didn't find any mention of pT."



"         MSTU(33)=1
          CALL PYROBO(N+1,N+NP+1,0D0,-PHI,0D0,0D0,0D0)
          THE=PYANGL(P(N+NP+1,3),P(N+NP+1,1))
          CALL PYROBO(N+1,N+NP+1,-THE,0D0,0D0,0D0,0D0)
        ENDIF
 
C...Find and order particles with highest p (pT for major).
        DO 110 ILF=N+NP+4,N+NP+MSTU(44)+4
          P(ILF,4)=0D0
  110   CONTINUE
        DO 160 I=N+1,N+NP
          IF(ILD.EQ.2) P(I,4)=SQRT(P(I,1)**2+P(I,2)**2)
          DO 130 ILF=N+NP+MSTU(44)+3,N+NP+4,-1
            IF(P(I,4).LE.P(ILF,4)) GOTO 140
"

The *reference* implementation is doing arrangement by pT, so it is the
expected behavior. Otherwise there will be a discrepancy between
calculated values.




2)
"
 what in particular seems suspicious 
> about the axis inversion"
That is the only case where the x is checked. In all other cases that is
z.


3) Also, same code looks very strange in other places:

"
vector<Vector3> p = momenta;
assert(p.size() >= 3);
unsigned int n = 3;
if (p.size() == 3) n = 3;
"

Andrii




On Tue, 2015-09-08 at 17:12 +0100, Andy Buckley wrote:
> On 08/09/15 16:20, Andrii Verbytskyi wrote:
> > Dear Andy,
> > 1) recently after looking into the thrust routine in Rivet I've fount
> > that objects are arranged by |p| and not p_T as in the referenced Pythia
> > code (line 77228 of pythia6428.f).
> > In some cases it produces unexpected outputs.
> > 2) The line
> > if (axis.x() < 0) axis = -axis;
> > in the same code looks very suspicious. Is that a bug or a feature?
> 
> Hi Andrii,
> 
> Thrust is principally an e+e- observable, so it makes sense that pT is 
> not prioritised. If you want transverse thrust then the input pz values 
> need to be set to 0 (or approximately zero -- I forget if the numerics 
> require some z component). I checked the header and documentation and 
> didn't find any mention of pT... did I miss something?
> 
> I've not looked into the code, but what in particular seems suspicious 
> about the axis inversion? This looks to me just like a way to 
> consistently orient the resulting set of principle axes. No physics 
> should depend on the absolute sign of the axis direction.
> 
> Andy
>