00001 // -*- C++ -*- 00002 #include "Rivet/Rivet.hh" 00003 #include "Rivet/Projections/DISLepton.hh" 00004 #include "Rivet/Cmp.hh" 00005 #include "Rivet/Tools/ParticleIdUtils.hh" 00006 00007 namespace Rivet { 00008 00009 int DISLepton::compare(const Projection& p) const { 00010 const DISLepton& other = pcast<DISLepton>(p); 00011 return 00012 mkNamedPCmp(other, "Beam") || 00013 mkNamedPCmp(other, "FS"); 00014 } 00015 00016 00017 void DISLepton::project(const Event& e) { 00018 const ParticlePair& inc = applyProjection<Beam>(e, "Beam").beams(); 00019 00020 Particle inLepton; 00021 00022 bool firstIsLepton = PID::isLepton(inc.first.pdgId()); 00023 bool secondIsLepton = PID::isLepton(inc.second.pdgId()); 00024 00025 if(firstIsLepton && !secondIsLepton){ 00026 _incoming = inc.first; 00027 }else if(!firstIsLepton && secondIsLepton){ 00028 _incoming = inc.second; 00029 }else{ 00030 //eek! 00031 throw Error("DISLepton projector could not find the correct beam. "); 00032 } 00033 00034 _sign = (_incoming.momentum().pz() > 0.0)? 1.0: -1.0; 00035 long id = _incoming.pdgId(); 00036 00037 double pzMax = -1000000000.0; 00038 00039 const FinalState& fs = applyProjection<FinalState>(e, "FS"); 00040 foreach (const Particle& p, fs.particles()) { 00041 double pz = _sign * p.momentum().pz(); 00042 if(p.pdgId() == id && pz > pzMax){ 00043 _outgoing = p; 00044 pzMax = pz; 00045 } 00046 } 00047 00048 if (!_outgoing.hasGenParticle()) { 00049 throw Error("DISLepton projector could not find the scattered lepton."); 00050 } 00051 } 00052 }