[Rivet] Bug in jet definition

[Andy Buckley]

Hi Roman,

Good: glad I understood. Thanks for bringing this up: we can definitely 
improve the interfaces and behaviours and should do so -- backward 
compatibility on things that probably shouldn't be used is a moot point ;)

For your interest, the reason for those methods on Rivet::Jet is that 
that CDF 2001 analysis was the first one that I wrote, and we needed a 
jet class. Fastjet didn't exist, and Rick's track-jet algorithm wasn't 
public, so I wrote a projection to implement it (or at least an 
interpretation of the ambiguous way that the paper describes it). I 
needed an output data structure, hence Rivet::Jet. But since we now use 
FastJet for everything (indeed, many of the plugins in FastJet came via 
the old Rivet implementations), some of those not-quite-duplicate 
methods should certainly be removed.

Cheers,
Andy


On 18/11/11 14:58, Roman Kogler wrote:
> Hi Andy,
>
> I looked at a handful events with the E-scheme and you are right, the
> mismatch between jets from Rivet and FastJet is gone. I see, it's the
> rescaling of particles to make them massless in FastJet which causes the
> difference. I think that Rivet should follow FastJet in this respect.
>
> I can understand that there is a good reason for having the Rivet::Jet
> interface. However, I'm not sure why you have to recalculate the jet
> four-vector after FastJet has already done this work for you. When I use
> jets, I expect that the four-vector of the jet is calculated according
> to the prescription of the recombination scheme which I pass to FastJet
> in the constructor:
> addProjection(FastJets(dfs, fastjet::kt_algorithm, fastjet::pt_scheme,
> 1.), "FastJet");
> This is exactly what happens in the data analysis, so it's best to use
> FastJet's calculation. Also I think it's safe for you to assume that any
> analysis that uses FastJet, will use its implementation of the
> recombination schemes. By the way, the recombination scheme has to be
> defined *before* the jet finding is performed, because it influences the
> jet finding itself.
>
> Concerning the pt-weighted jet properties in Rivet, I can't help you.
> I'm not sure where these are used and why they have been implemented.
> It's of course always useful to have additional methods available, but
> to be honest, I find the methods from Rivet::Jet a bit confusing. For
> example, what's the difference between
> Jet::ptWeightedEta() and Jet::eta()? When using the pt-recombination
> scheme there should not be a difference. However, we found out that
> there is a difference because of a different treatment of the particle
> masses, this can be confusing. Anyways, maybe you would still like to
> keep both methods for backward compatibility. I saw that there is only
> one analysis which uses the Jet::ptSum() method: CDF_2001_S4751469
> <http://projects.hepforge.org/rivet/code/dev/a00315_source.html#l00086>.
> Do you know which recombination scheme this analysis exactly used?
>
> In any event, I think the important point is that the jet four-vector
> should be filled by FastJet.
>
> Cheers,
> Roman
>
>
>
>
>
>
> On 11/18/11 0:40, Andy Buckley wrote:
>> Hi Roman,
>>
>> I just had a bit of a look at this -- to check that I've understood
>> this correctly, is this mismatch issue only seen when using a jet
>> recombination scheme other than E (i.e. pure 4-vector addition)?
>>
>> Cross-checking FastJet's manual re. the pT scheme and Rivet's
>> pT-weighted jet properties, there seems to be a difference in that the
>> FastJet one rescales momenta to be massless before combining. Rivet
>> doesn't do that: we just do this
>>
>> foreach (const FourMomentum& p, momenta()) {
>> double pt = p.pT();
>> ptsum += pt;
>> ptwetasum += pt * p.pseudorapidity();
>> ptwdphisum += pt * mapAngleMPiToPi(phi0 - p.azimuthalAngle());
>> }
>> _ptWeightedEta = ptwetasum/ptsum;
>> [...]
>>
>> Maybe that's definitely wrong -- if so we can fix it. Or just remove
>> it / base it on the accumulated 4-vector. There is probably work that
>> can be done to improve the Rivet::Jet interface, which was not very
>> carefully designed and doesn't provide options for recombination other
>> than the E scheme (but does have other reasons to exist).
>>
>> Anyway, if I've understood correctly this is an area where I'd like to
>> hear suggestions about how we can do better, but I'm breathing a sign
>> of relief because I think all the "official" analyses that I know are
>> using E-recombination and should be unaffected.
>>
>> Andy
>>
>>
>> On 14/11/11 14:12, Roman Kogler wrote:
>>> Dear Rivet developers,
>>>
>>> I believe I found a bug in the jet definition of Rivet. It's a quite
>>> annoying feature and it took me quite some time to figure it out.
>>>
>>> First, I couldn't reproduce the results that come out of my jet analysis
>>> in H1. To find the difference I put some Sherpa events into the format
>>> we use in H1 for our analyses. Then I started to compare the output
>>> even-by-event, using the same events in my H1 analysis and the Rivet
>>> analysis I wrote. The boost of the HFS to the Breit frame is identical,
>>> as well as the number of jets found. Even the jet constituents of each
>>> jet are the same in both analyses. So finally I tracked the difference
>>> to the jet four-vector. It seems that Rivet overwrites the jet
>>> four-vector as calculated by FastJet. To clarify that I looked into the
>>> FastJets projection, Projections/FastJets.C. I put some output
>>> statements into the method _pseudojetsToJets, the output is labelled
>>> "from FastJet" below (all code can be found at the end of this email).
>>> Also, I put some output into my analysis using the standard
>>> jet::momentum() method ("from Rivet"). I then recalculate the jet
>>> four-vector with the pt-recombination scheme ("my calculation"). Here
>>> are the results for one event:
>>>
>>> jet 0: from FastJet, Pt = 14.6163 phi = 6.27012 eta = 0.406647
>>> jet 1: from FastJet, Pt = 14.3426 phi = 3.08307 eta = 0.240164
>>>
>>> jet 0: from Rivet, Pt = 14.0764 phi = 6.27134 eta = 0.443595, pt
>>> weighted phi = -0.0106254 pt weighted eta = 0.406647
>>> jet 1: from Rivet, Pt = 13.5777 phi = 3.08787 eta = 0.259143, pt
>>> weighted phi = 3.09267 pt weighted eta = 0.240164
>>>
>>> jet 0: my calculation: pt = 14.6163 phi = 6.27012 eta = 0.406647
>>> jet 1: my calculation: pt = 14.3426 phi = 3.08307 eta = 0.240164
>>>
>>> The jet four-vectors from FastJet agree with the ones using my
>>> implementation of the pt-recombination scheme. Also, these numbers are
>>> identical with the numbers I get using the H1 analysis framework (with
>>> its own interface to FastJet). But the jets stored by Rivet have
>>> different pt! How could this happen? It is really quite a dangerous bug,
>>> I believe. Please have a look at it soon.
>>>
>>> Cheers,
>>> Roman
>>>
>>>
>>> P.S. Here is the code I used for the output above:
>>>
>>>
>>> -----------------------------------------------
>>> Projections/FastJets.C:
>>>
>>> Jets FastJets::_pseudojetsToJets(const PseudoJets& pjets) const {
>>> Jets rtn;
>>> int i=0;
>>> foreach (const fastjet::PseudoJet& pj, pjets) {
>>> cout << "jet " << i++ << ": ";
>>> cout << "from FastJet, Pt = " << pj.perp() << " phi = " << pj.phi() << "
>>> eta = " << pj.rap() << endl;
>>> ...
>>> }
>>>
>>>
>>> -----------------------------------------------
>>> H1 jet analysis:
>>>
>>> void init() {
>>> ...
>>> addProjection(FastJets(dfs, fastjet::kt_algorithm, fastjet::pt_scheme,
>>> 1.), "FastJet");
>>> ...
>>> }
>>>
>>> void analyze(const Event& event) {
>>>
>>> const FastJets & fj = applyProjection<FastJets>(event, "FastJet");
>>> Jets jets = fj.jetsByPt(3., 80.);
>>>
>>> vector<FourMomentum> jet_vecs;
>>> int njet = 0;
>>>
>>> foreach (const Jet& j, jets){
>>>
>>> cout << "jet " << njet++ << ": ";
>>> cout << "from Rivet, Pt = " << j.momentum().pT() << " phi = " <<
>>> j.momentum().phi() << " eta = " << j.momentum().eta()
>>> << " pt weighted phi = " << j.ptWeightedPhi() << " pt weighted eta = "
>>> << j.ptWeightedEta() << endl;
>>>
>>> // recalculate the jet four-vector with my own implementation of the
>>> pt-recombination scheme
>>> const vector<Particle> parts = j.particles();
>>> FourMomentum jet_vec(0, 0, 0, 0);
>>>
>>> foreach(const Particle& jp, parts){
>>>
>>> // scale particle momtenta to make them massless
>>> double pp = sqrt(jp.momentum().px()*jp.momentum().px() +
>>> jp.momentum().py()*jp.momentum().py() +
>>> jp.momentum().pz()*jp.momentum().pz() );
>>> FourMomentum pv(pp, jp.momentum().px(), jp.momentum().py(),
>>> jp.momentum().pz());
>>> if (ij==0){
>>> jet_vec = pv;
>>> } else {
>>> jet_vec = MergeParts(pv, jet_vec);
>>> }
>>> ++ij;
>>> }
>>>
>>> cout << "my calculation: pt = " << jet_vec.pT() << " phi = " <<
>>> jet_vec.phi() << " eta = " << jet_vec.eta() << endl;
>>> }
>>>
>>>
>>> FourMomentum MergeParts(FourMomentum a, FourMomentum b)
>>> {
>>> // merge momenta a and b with pt-recombination scheme
>>>
>>> double pt = a.pT() + b.pT();
>>> double phi_a = a.phi();
>>> double phi_b = b.phi();
>>> if (phi_a - phi_b > fastjet::pi) phi_b += fastjet::twopi;
>>> if (phi_a - phi_b < -fastjet::pi) phi_b -= fastjet::twopi;
>>> double phi_sum = a.pT()*phi_a + b.pT()*phi_b;
>>> double phi = phi_sum / pt;
>>> double eta_sum = a.pT()*a.pseudorapidity() + b.pT()*b.pseudorapidity();
>>> double eta = eta_sum / pt;
>>>
>>> FourMomentum v(pt*cosh(eta), pt*cos(phi), pt*sin(phi), pt*sinh(eta));
>>> return v;
>>>
>>> }
>>>
>>>
>>
>>


-- 
Dr Andy Buckley
SUPA Advanced Research Fellow
Particle Physics Experiment Group, University of Edinburgh

The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.