[Rivet] Questions on Rivet for ep DIS

[Roman Kogler]

Hi Hendrik et al.,

the list from Thomas is a follow-up on a discussion I had some time ago 
with Andy. At that time, I started to implement a DIS high-Q2 jet 
analysis in Rivet and found some show-stoppers:

1) There is a bug in DISLepton.cc. The electron (or positron) with the 
maximum p_z (measured in the direction of the electron beam) is taken to 
be the scattered electron (positron) in the event. While this holds 
approximately at low Q2, it's completely wrong at high Q2, where the 
incoming electron can be back-scattered and actually have opposite p_z 
w.r.t. the electron beam. We had a lengthy discussion about how the 
scattered electron should be defined, since there is always an ambiguity 
if there are multiple electrons in the event. However, the H1 and ZEUS 
experiments always correct for a possible identification inefficiency of 
the scattered electron, so the scattered electron has to be found with a 
100% efficiency in Rivet. I just checked on hepforge and I think this 
bug is still present in the current version.

2) There is a bug in the jet definition if the pt-recombination scheme 
is used. Here Rivet overwrites the jet fourvector of FastJet and 
recomputes the fourvector in a wrong way. Below are some details. I'm 
not sure if this has been fixed.

3) I needed the hadronic final state in the Breit frame of reference for 
the jet clustering. This can be obtained by using the DISFinalState 
projection and having an additional argument in the constructor, where 
one can choose between the HCM and Breit frame. Having a look at the 
code documentation, this seems to be in now.

Could you please tell us what's the status of points 1 and 2?

Cheers,
   Roman




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.




On 6/11/12 15:26, Andy Buckley wrote:
> Hi Thomas,
>
> Just replying to include the Rivet developer mailing list in the
> discussion (hopefully one of the other developers will be able to answer
> properly sooner than I can)
>
> Andy
>
>
> On 11/06/12 08:57, Thomas Schoerner-Sadenius wrote:
>> Dear Andy, (ss Roman, Simon),
>>
>> in a recent discussion with Roman on new Rivet routines for HERA jet
>> physics Roman mentioned a few problems that existed some time ago, and
>> he was unsure about the current status. He suggested to contact you - so
>> maybe you can help with the following issues? Thanks a lot in advance!
>>
>> - Electron finding for very high Q2 values. According to Roman the issue
>> here was that the electron was picked according to the maximum pz -
>> which fails for high Q2. Has that been fixed ?
>> - A problem in the merging of particles for the pT reco scheme in the kT
>> finder?
>> - And … no, I think that was it. Roman, did I forget anything?
>>
>> I am very much looking forward to hearing from you!
>>
>> Kind regards,
>>
>> Thomas
>>
>>
>>
>> -- 
>>
>> ============================
>>
>> Thomas Schoerner-Sadenius
>> Analysis Centre, Helmholtz Alliance ‘Physics at the Terascale’
>> DESY/FH-CMS, Notkestr. 85, D-22607 Hamburg, Germany
>> e-mail: thomas.schoerner at desy.de<mailto:thomas.schoerner at desy.de>
>> phone: +49 40 8998 3429        fax: +49 40 8998 4304
>> http://www.desy.de/~schorner
>>
>> http://www.terascale.de
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.hepforge.org/lists-archive/rivet/attachments/20120612/16ed8eed/attachment.html>