[Rivet-svn] r1679 - in trunk: include/Rivet/Math include/Rivet/Projections src/Projections

[blackhole at projects.hepforge.org]

Author: buckley
Date: Mon Jul 13 21:00:03 2009
New Revision: 1679

Log:
Fixing definition of HCM Lorentz Transform to put the virtual photon along the +z direction for the HCM frame, and the -z direction for the Breit frame

Modified:
   trunk/include/Rivet/Math/Vector3.hh
   trunk/include/Rivet/Math/Vector4.hh
   trunk/include/Rivet/Projections/DISKinematics.hh
   trunk/src/Projections/DISKinematics.cc
   trunk/src/Projections/FinalStateHCM.cc
   trunk/src/Projections/TotalVisibleMomentum.cc

Modified: trunk/include/Rivet/Math/Vector3.hh
==============================================================================
--- trunk/include/Rivet/Math/Vector3.hh	Mon Jul 13 20:58:42 2009	(r1678)
+++ trunk/include/Rivet/Math/Vector3.hh	Mon Jul 13 21:00:03 2009	(r1679)
@@ -79,11 +79,8 @@
     }
 
     double angle(const Vector3& v) const {
-      
-      double localDotOther = unit().dot(v.unit());
-      
+      double localDotOther = unit().dot(v.unit());      
       if(Rivet::isZero(localDotOther - 1.0)) return 0.0;
-      
       return acos( localDotOther );
     }
 

Modified: trunk/include/Rivet/Math/Vector4.hh
==============================================================================
--- trunk/include/Rivet/Math/Vector4.hh	Mon Jul 13 20:58:42 2009	(r1678)
+++ trunk/include/Rivet/Math/Vector4.hh	Mon Jul 13 21:00:03 2009	(r1679)
@@ -27,20 +27,20 @@
 
     template<typename V4>
     FourVector(const V4& other) {
-      this->t(other.t());
-      this->x(other.x());
-      this->y(other.y());
-      this->z(other.z());
+      this->setT(other.t());
+      this->setX(other.x());
+      this->setY(other.y());
+      this->setZ(other.z());
     }
 
     FourVector(const Vector<4>& other) 
     : Vector<4>(other) { }
 
     FourVector(const double t, const double x, const double y, const double z) {
-      this->t(t);
-      this->x(x);
-      this->y(y);
-      this->z(z);
+      this->setT(t);
+      this->setX(x);
+      this->setY(y);
+      this->setZ(z);
     }
 
     virtual ~FourVector() { }
@@ -50,10 +50,10 @@
     const double x() const { return get(1); }
     const double y() const { return get(2); }
     const double z() const { return get(3); }
-    FourVector& t(const double t) { set(0, t); return *this; }
-    FourVector& x(const double x) { set(1, x); return *this; }
-    FourVector& y(const double y) { set(2, y); return *this; }
-    FourVector& z(const double z) { set(3, z); return *this; }
+    FourVector& setT(const double t) { set(0, t); return *this; }
+    FourVector& setX(const double x) { set(1, x); return *this; }
+    FourVector& setY(const double y) { set(2, y); return *this; }
+    FourVector& setZ(const double z) { set(3, z); return *this; }
 
     double invariant() const {
       return t()*t() - x()*x() - y()*y() - z()*z();
@@ -311,20 +311,20 @@
 
     template<typename V4>
     FourMomentum(const V4& other) {
-      this->E(other.t());
-      this->px(other.x());
-      this->py(other.y());
-      this->pz(other.z());
+      this->setE(other.t());
+      this->setPx(other.x());
+      this->setPy(other.y());
+      this->setPz(other.z());
     }
 
     FourMomentum(const Vector<4>& other) 
       : FourVector(other) { }
 
     FourMomentum(const double E, const double px, const double py, const double pz) {
-      this->E(E);
-      this->px(px);
-      this->py(py);
-      this->pz(pz);
+      this->setE(E);
+      this->setPx(px);
+      this->setPy(py);
+      this->setPz(pz);
     }
 
     ~FourMomentum() {}
@@ -346,16 +346,16 @@
     double pz() const { return z(); }
 
     /// Set energy \f$ E \f$ (time component of momentum).
-    FourMomentum& E(double E)   { t(E); return *this; }
+    FourMomentum& setE(double E)   { setT(E); return *this; }
 
     /// Set x-component of momentum \f$ p_x \f$.
-    FourMomentum& px(double px) { x(px); return *this; }
+    FourMomentum& setPx(double px) { setX(px); return *this; }
 
     /// Set y-component of momentum \f$ p_y \f$.
-    FourMomentum& py(double py) { y(py); return *this; }
+    FourMomentum& setPy(double py) { setY(py); return *this; }
 
     /// Set z-component of momentum \f$ p_z \f$.
-    FourMomentum& pz(double pz) { z(pz); return *this; }
+    FourMomentum& setPz(double pz) { setZ(pz); return *this; }
 
     /// Get squared mass \f$ m^2 = E^2 - p^2 \f$ (the Lorentz self-invariant).
     double mass2() const { 

Modified: trunk/include/Rivet/Projections/DISKinematics.hh
==============================================================================
--- trunk/include/Rivet/Projections/DISKinematics.hh	Mon Jul 13 20:58:42 2009	(r1678)
+++ trunk/include/Rivet/Projections/DISKinematics.hh	Mon Jul 13 21:00:03 2009	(r1679)
@@ -19,11 +19,7 @@
         
     /// The default constructor.
     DISKinematics() 
-    :  _theQ2(-1.0), 
-       _theW2(-1.0), 
-       _theX(-1.0), 
-       _theY(-1.0), 
-       _theS(-1.0) 
+      : _theQ2(-1.0), _theW2(-1.0), _theX(-1.0), _theY(-1.0), _theS(-1.0) 
     {
       setName("DISKinematics");
       //addBeamPair(ANY, hadid);

Modified: trunk/src/Projections/DISKinematics.cc
==============================================================================
--- trunk/src/Projections/DISKinematics.cc	Mon Jul 13 20:58:42 2009	(r1678)
+++ trunk/src/Projections/DISKinematics.cc	Mon Jul 13 21:00:03 2009	(r1679)
@@ -10,12 +10,10 @@
 
 
   void DISKinematics::project(const Event& e) {
-    const DISLepton& dislep = applyProjection<DISLepton>(e, "Lepton");
+    // Identify beam hadron
     const ParticlePair& inc = applyProjection<Beam>(e, "Beam").beams();
-
     bool firstIsHadron  = PID::isHadron(inc.first.pdgId());
-    bool secondIsHadron = PID::isHadron(inc.second.pdgId());
-    
+    bool secondIsHadron = PID::isHadron(inc.second.pdgId());    
     if (firstIsHadron && !secondIsHadron) {
       _inHadron = inc.first;
     } else if (!firstIsHadron && secondIsHadron) {
@@ -25,6 +23,8 @@
       throw Error("DISKinematics projector could not find the correct beam hadron");
     }
 
+    // Get the DIS lepton and store some of its properties
+    const DISLepton& dislep = applyProjection<DISLepton>(e, "Lepton");
     const FourMomentum pLepIn = dislep.in().momentum();
     const FourMomentum pLepOut = dislep.out().momentum();
     const FourMomentum pHad = _inHadron.momentum();
@@ -37,34 +37,36 @@
     _theS = invariant(pLepIn + pHad);
 
     // Calculate boost vector for boost into HCM-system
-    _hcm.setBoost(-tothad.boostVector());
+    LorentzTransform tmp;
+    tmp.setBoost(-tothad.boostVector());
 
-    // Boost the gamma and lepton
-    FourMomentum pGammaHCM = _hcm.transform(pGamma);
-
-    // Rotate so the photon is in x-z plane
-    _hcm.preMult(Matrix3(Vector3::mkZ(), -pGammaHCM.azimuthalAngle()));
-    pGammaHCM = _hcm.transform(pGamma);
+    // Rotate so the photon is in x-z plane in HCM rest frame
+    FourMomentum pGammaHCM = tmp.transform(pGamma);
+    tmp.preMult(Matrix3(Vector3::mkZ(), -pGammaHCM.azimuthalAngle()));
+    pGammaHCM = tmp.transform(pGamma);
     assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkY())));
 
-    // Rotate so the photon is along the negative z-axis
-    _hcm.preMult(Matrix3(Vector3::mkY(), pi - pGammaHCM.polarAngle()));
-
-    // Check that final HCM photon lies along -ve z as expected
-    pGammaHCM = _hcm.transform(pGamma);
-    assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkX())) &&
-           isZero(dot(pGammaHCM.vector3(), Vector3::mkY())));
-    assert(isZero(angle(pGammaHCM.vector3(), -Vector3::mkZ())));
+    // Rotate so the photon is along the positive z-axis
+    const double rot_angle = pGammaHCM.polarAngle() * (pGammaHCM.px() >= 0 ? -1 : 1);
+    tmp.preMult(Matrix3(Vector3::mkY(), rot_angle));
+    // Check that final HCM photon lies along +ve z as expected
+    pGammaHCM = tmp.transform(pGamma);
+    assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkX())));
+    assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkY())));
+    assert(isZero(angle(pGammaHCM.vector3(), Vector3::mkZ())));
 
     // Finally rotate so outgoing lepton at phi = 0
-    FourMomentum pLepOutHCM = _hcm.transform(pLepOut);
-    _hcm.preMult(Matrix3(Vector3::mkZ(), -pLepOutHCM.azimuthalAngle()));
-    pLepOutHCM = _hcm.transform(pLepOut);
-    assert(isZero(pLepOutHCM.azimuthalAngle()));
+    FourMomentum pLepOutHCM = tmp.transform(pLepOut);
+    tmp.preMult(Matrix3(Vector3::mkZ(), -pLepOutHCM.azimuthalAngle()));
+    assert(isZero(tmp.transform(pLepOut).azimuthalAngle()));
+    _hcm = tmp;
 
-    // Boost to Breit frame    
+    // Boost to Breit frame (use opposite convention for photon --- along *minus* z)
+    tmp.preMult(Matrix3(Vector3::mkX(), PI));
     const double bz = 1 - 2*x();
-    _breit = LorentzTransform(Vector3::mkZ() * bz).combine(_hcm);
+    _breit = LorentzTransform(Vector3::mkZ() * bz).combine(tmp);
+    assert(isZero(angle(_breit.transform(pGamma).vector3(), -Vector3::mkZ())));
+    assert(isZero(_breit.transform(pLepOut).azimuthalAngle()));
   }
 
 

Modified: trunk/src/Projections/FinalStateHCM.cc
==============================================================================
--- trunk/src/Projections/FinalStateHCM.cc	Mon Jul 13 20:58:42 2009	(r1678)
+++ trunk/src/Projections/FinalStateHCM.cc	Mon Jul 13 21:00:03 2009	(r1679)
@@ -5,9 +5,9 @@
 
 namespace Rivet {
 
+
   int FinalStateHCM::compare(const Projection& p) const {
-    return \
-      mkNamedPCmp(p, "Kinematics");
+    return mkNamedPCmp(p, "Kinematics");
   }
 
 
@@ -33,4 +33,5 @@
     }
   }
 
+
 }

Modified: trunk/src/Projections/TotalVisibleMomentum.cc
==============================================================================
--- trunk/src/Projections/TotalVisibleMomentum.cc	Mon Jul 13 20:58:42 2009	(r1678)
+++ trunk/src/Projections/TotalVisibleMomentum.cc	Mon Jul 13 21:00:03 2009	(r1679)
@@ -3,9 +3,9 @@
 #include "Rivet/Projections/TotalVisibleMomentum.hh"
 #include "Rivet/Cmp.hh"
 
-
 namespace Rivet {
 
+
   int TotalVisibleMomentum::compare(const Projection& p) const {
     return mkNamedPCmp(p, "FS");
   }
@@ -13,17 +13,18 @@
 
   void TotalVisibleMomentum::project(const Event& e) {
     _momentum = FourMomentum();
-    _momentum.px(0.0).py(0.0).pz(0.0).E(0.0);
     _set = 0.0;
 
     // Project into final state
     const FinalState& fs = applyProjection<FinalState>(e, "FS");
 
     // Get hadron and charge info for each particle, and fill counters appropriately
-    for (ParticleVector::const_iterator p = fs.particles().begin(); p != fs.particles().end(); ++p) {
-      _momentum += p->momentum();
-      _set += pT(p->momentum());
+    foreach (const Particle& p, fs.particles()) {
+      const FourMomentum& mom = p.momentum();
+      _momentum += mom;
+      _set += mom.pT();
     }
   }
 
+
 }