[Rivet] Bug in jet definition

[Roman Kogler]

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;

   }