[Rivet] Jets in DIS

[Roman Kogler]

Dear Andy,

I checked out the new Rivet version, thanks for including the boost of 
the final state to the Breit frame. I'm happy to see that working.

However, there is a bug in Rivet when looking for the scattered electron 
in DIS. The algorithm (DISLepton.cc) looks for the electron (or 
positron) in the event with the highest value of pz. This works at low 
Q2 rather well, but at high Q2 this algorithm fails. Here is why:
At HERA the +pz direction is defined to be the proton beam direction. So 
the electron beam has pz=-27.6. At low Q2, where the electron is 
scattered under small angles, the 'signed value' of the scattered 
electron is then +pz'. However, at high Q2 the electron can be scattered 
perpendicular to the beam (pz'=0) or in the forward direction (-pz'). So 
if there is any other electron in the event present, which goes in the 
-pz direction, it will be misidentified as the scattered electron. Below 
is an example for an event generated with Q2~17000, where the scattered 
electron has been misidentified and the calculated value of Q2 in Rivet 
is 72.8.

Cheers,
    Roman


DISLepton: searching for ID -11 with pz > -1e+09
part 0 id = -11 px = -55.7936 py = 67.9235 pz = 73.3469 E = 114.483
part 0 id = -11 signed_pz = -73.3469
   ---> taking 0 as outgoing lepton
part 1 id = 211 px = 4.3922 py = -5.36174 pz = 21.1347 E = 22.2426
part 2 id = 2212 px = 19.1258 py = -22.3493 pz = 85.5571 E = 90.4775
part 3 id = -2212 px = 7.27322 py = -8.24468 pz = 32.7272 E = 34.5373
part 4 id = 2112 px = -0.258489 py = 0.0192201 pz = 140.113 E = 140.116
part 5 id = 22 px = 0.0752175 py = -0.0810541 pz = 0.527222 E = 0.538693
part 6 id = -211 px = 2.95982 py = -3.164 pz = 14.1068 E = 14.7578
part 7 id = 130 px = 0.251775 py = -0.163703 pz = 183.251 E = 183.252
part 8 id = 211 px = 6.86782 py = -8.30335 pz = 33.2207 E = 34.9249
part 9 id = -211 px = 2.33043 py = -2.5695 pz = 10.0595 E = 10.6417
part 10 id = 22 px = 0.693957 py = -0.670989 pz = 3.14172 E = 3.28667
part 11 id = 22 px = 1.11994 py = -1.31727 pz = 5.41137 E = 5.68088
part 12 id = 130 px = -0.151476 py = -0.389603 pz = 18.067 E = 18.0787
part 13 id = -321 px = -0.516672 py = -0.553757 pz = 19.2324 E = 19.2536
part 14 id = 211 px = 0.0829556 py = 0.0993597 pz = 4.24924 E = 4.2535
part 15 id = 130 px = 0.375512 py = -0.621349 pz = 14.7582 E = 14.7844
part 16 id = 22 px = -0.0301481 py = -0.00922794 pz = 0.33583 E = 0.337307
part 17 id = 22 px = 0.340224 py = -0.445809 pz = 1.41281 E = 1.52004
part 18 id = 11 px = 0.075713 py = -0.141711 pz = 0.556192 E = 0.578934
part 19 id = -11 px = 0.124663 py = -0.154993 pz = 0.645015 E = 0.674987
part 19 id = -11 signed_pz = -0.645015
   ---> taking 19 as outgoing lepton
part 20 id = 22 px = 0.04499 py = -0.12844 pz = 0.397113 E = 0.419785
part 21 id = 22 px = 0.827752 py = -1.19367 pz = 4.91671 E = 5.1268
part 22 id = 22 px = 0.748384 py = -1.12375 pz = 4.89815 E = 5.08082
part 23 id = 22 px = -0.00267821 py = -0.00925982 pz = 0.0143985 E = 
0.0173273
part 24 id = 130 px = 0.287098 py = 0.19742 pz = 17.3042 E = 17.3149
part 25 id = 130 px = 0.122298 py = -0.0634414 pz = 25.2863 E = 25.2916
part 26 id = 22 px = 0.0347235 py = 0.144282 pz = 1.65602 E = 1.66266
part 27 id = 22 px = -0.0248541 py = 0.146226 pz = 3.67157 E = 3.67456
part 28 id = 211 px = 0.954048 py = -0.963005 pz = 3.13684 E = 3.42006
part 29 id = 211 px = 0.235199 py = -0.343671 pz = 1.4502 E = 1.51525
part 30 id = -12 px = 1.57995 py = -1.37946 pz = 7.29499 E = 7.59052
part 31 id = 11 px = 0.471629 py = -0.856397 pz = 2.71535 E = 2.88599
part 32 id = 211 px = -0.0294085 py = -0.547087 pz = 0.758253 E = 0.94583
part 33 id = -211 px = 0.0861909 py = -0.0196154 pz = 0.452828 E = 0.482023
part 34 id = 211 px = 0.225543 py = -0.161993 pz = 0.915438 E = 0.966756
part 35 id = -211 px = 0.571015 py = -0.738361 pz = 5.30456 E = 5.38786
part 36 id = 211 px = 0.129393 py = -0.267621 pz = 1.02406 E = 1.07542
part 37 id = -211 px = 1.77415 py = -1.73558 pz = 7.05777 E = 7.48274
part 38 id = 22 px = 0.205595 py = -0.0964995 pz = 1.00644 E = 1.03175
part 39 id = 22 px = -0.00465838 py = -0.055959 pz = 0.356084 E = 0.360484
part 40 id = 22 px = 0.567172 py = -0.788086 pz = 3.53189 E = 3.66293
part 41 id = 22 px = 0.406857 py = -0.734593 pz = 2.77624 E = 2.90046
part 42 id = 211 px = 0.374462 py = -0.325563 pz = 2.74971 E = 2.7976
part 43 id = -211 px = 0.206282 py = -0.339635 pz = 0.964825 E = 1.05274
part 44 id = -211 px = 1.22403 py = -0.939285 pz = 24.7002 E = 24.7488
part 45 id = 211 px = -0.401421 py = 0.084246 pz = 17.7685 E = 17.7737
part 46 id = 22 px = 0.657991 py = -0.669357 pz = 29.0299 E = 29.0451
part 47 id = 22 px = 0.30748 py = 0.0117005 pz = 5.16715 E = 5.17631
part 48 id = 22 px = -0.0346924 py = -0.117068 pz = 6.57012 E = 6.57125
part 49 id = 321 px = -0.74251 py = -0.217304 pz = 43.1924 E = 43.2021
part 50 id = -321 px = -0.140918 py = -0.269146 pz = 4.47808 E = 4.51544
fulfilled phase space cuts, Q2 = 72.8641 y = 0.999457 weight = 51.0074
scat. elec: E = 0.674987 theta = 0.299122

Here is the code that I used the generate this output. It's the code of 
DISLepton.cc with some output statements added.

     cout << "DISLepton: searching for ID " << id << " with pz > " << 
pzMax << endl;
     const FinalState& fs = applyProjection<FinalState>(e, "FS");
     int ipart = 0;
     foreach (const Particle& p, fs.particles()) {

       cout << "part " << ipart << " id = " << p.pdgId()
<< " px = " << p.momentum().px()
<< " py = " << p.momentum().py()
<< " pz = " << p.momentum().pz()
<< " E = " << p.momentum().E() << endl;
       double pz = _sign * p.momentum().pz();
       if(p.pdgId() == id && pz > pzMax){
         cout << "part " << ipart << " id = " << p.pdgId() << " 
signed_pz = " << _sign * p.momentum().pz() << endl;
         cout << "  ---> taking " << ipart << " as outgoing lepton " << 
endl;
         _outgoing = p;
         pzMax = pz;
       }
       ++ipart;
     }




On 8/15/11 10:17, Andy Buckley wrote:
> Hi Roman,
>
> I'm no DIS expert, but I think that the FinalStateHCM should do 
> something like that -- is it not acceptable?
>
> If that is not appropriate, you don't absolutely need to use the 
> "central registry" aspect of the projection system: you can make your 
> list of particles, boosted and filtered however you want, in the 
> analyze() method, and then use the calc() method of a local JetAlg 
> object to make the jets. Let us know if you need more help with this 
> direction: we've only recently started providing this sort of access 
> to projection features without having to register the projections -- 
> sometimes it's necessary, but the normal docs don't describe it!
>
> Cheers,
> Andy
>
>
> On 15/08/11 18:08, Roman Kogler wrote:
>> Dear Rivet authors and developers,
>>
>> I'm trying to implement a new jet analysis in Rivet. The data come from
>> H1. The analysis is not published yet, but we are close to the
>> publication. It's a measurement of inclusive jet, dijet and trijet cross
>> sections in DIS at high Q^2 (for more information have a look at
>> H1prelim-11-032).
>>
>> In this jet analysis jets are found in the Breit frame and I could not
>> find a way in Rivet how to do this. What I need is to boost the full
>> event to the Breit frame before the jet finding is performed. FastJet
>> should then be applied to all final state particles in the Breit frame.
>>
>> Is there an elegant solution to do this in Rivet?
>>
>> Thanks,
>> Roman
>>
>>
>
>