[Rivet] thrust

[Andy Buckley]

On 08/09/15 17:43, Andrii Verbytskyi wrote:
> 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.

First, that comment is in the Thrust.cc implementation file and not 
promised by the user-facing comments on the Thrust projection interface. 
Further, I've checked and in fact the comment is correct in that it 
implements the algorithm exactly as described in the Pythia manual, 
which specifically describes the ordering as being in absolute momentum. 
We don't promise to implement exactly the same details as what is 
implemented in Pythia, so unless you have proof that the Pythia manual 
description is incorrect then I think we should leave the Rivet 
implementation as it is. From a physics point of view it makes no sense 
to privilege the beam direction in e+e-.

> 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.

Again, we have not declared (and certainly not publicly) that Pythia's 
implementation is any sort of canonical reference for ours. From the 
comment in that code, it seems clear to me that absolute p is used in 
general, and probably that "pT" refers to the momentum transverse to the 
already-identified thrust axis (which makes sense for the thrust major). 
I would have to go through the code with a fine-tooth comb to reverse 
engineer exactly what it is doing as compared to the definition in the 
Pythia manual.

> 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.

It doesn't matter: we can orient the axes any way we like as long as we 
do it consistently. No axis is special in e+e-.

> 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;
> "

Agreed that it looks strange... in fact I think I remember this bit 
arising from step-wise removal of some earlier bad logic. But unless 
there's actually evidence of a physics bug I don't think there is any 
pressing need to change it.

Andy



> 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
>>
>
>


-- 
Dr Andy Buckley, Lecturer / Royal Society University Research Fellow
Particle Physics Expt Group, University of Glasgow