use covariant formulas for PhiR and Theta

Author: c-dilks

src/Dihadrons.cxx

@@ -5,8 +5,8 @@ ClassImp(Dihadron)
 
 Dihadron::Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K)
 {
-  const TLorentzVector &vecEB  = K->vecEleBeam;
-  const TLorentzVector &vecIB  = K->vecIonBeam;
+  const TLorentzVector &vecE   = K->vecEleBeam;
+  const TLorentzVector &vecP   = K->vecIonBeam;
   const TLorentzVector &vecL   = K->vecElectron;
   const TLorentzVector &vecQ   = K->vecQ;
   const TLorentzVector &vecW   = K->vecW;
@@ -23,42 +23,44 @@ Dihadron::Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K
   for(int h=0; h<2; h++)
     K->BoostToIonFrame(vecH[h],IvecH[h]);
 
-  // invariant mass, missing mass, Z, PhPerp
+  // invariant mass, missing mass, Z, PhPerp, Zeta
   Mh = vecPh.M();
   MX = TMath::Abs((vecW-vecPh).M());
-  Z  = vecIB.Dot(vecPh) / vecIB.Dot(vecQ);
+  Z  = vecP.Dot(vecPh) / vecP.Dot(vecQ);
   for(int h=0; h<2; h++)
-    hadZ[h] = vecIB.Dot(vecH[h]) / vecIB.Dot(vecQ);
+    hadZ[h] = vecP.Dot(vecH[h]) / vecP.Dot(vecQ);
   PhPerp = Kinematics::Reject(IvecPh.Vect(),IvecQ.Vect()).Mag();
-
-  // Theta
-  auto boost_com = -1 * vecPh.BoostVector(); // boost to dihadron CoM frame
-  TLorentzVector vecH_com[2];
-  for(int h=0; h<2; h++) {
-    vecH_com[h] = vecH[h];
-    vecH_com[h].Boost(boost_com);
-  }
-  Theta = Kinematics::AngleSubtend(vecH_com[0].Vect(),vecPh.Vect());
+  Zeta = 2 * vecR.Dot(vecP) / vecPh.Dot(vecP);
 
   // PhiH
   PhiH = Kinematics::PlaneAngle(
-      vecQ.Vect(), vecL.Vect(),
-      vecQ.Vect(), vecPh.Vect()
+      IvecQ.Vect(), IvecL.Vect(),
+      IvecQ.Vect(), IvecPh.Vect()
       );
 
-  // PhiR (calculated in ion rest frame)
-  TVector3 momHad_perp[2]; // perp-frame hadron 3-momenta
-  for(int h=0; h<2; h++)
-    momHad_perp[h] = Kinematics::Reject(IvecH[h].Vect(), IvecQ.Vect());
-  auto momRT = 1 / ( hadZ[0] + hadZ[1] ) *  
-    ( hadZ[1] * momHad_perp[0]  -  hadZ[0] * momHad_perp[1] );
+  // PhiR
+  Double_t coeff =
+    K->x * (Zeta*Mh*Mh - (vecH[0].M2() - vecH[1].M2()) ) /
+    ( K->Q2 * Z );
+  TLorentzVector vecRT = vecR - (Zeta/2)*vecPh + coeff*vecP;
+  TLorentzVector IvecRT;
+  K->BoostToIonFrame(vecRT,IvecRT);
   PhiR = Kinematics::PlaneAngle(
       IvecQ.Vect(), IvecL.Vect(),
-      IvecQ.Vect(), momRT
+      IvecQ.Vect(), IvecRT.Vect()
       );
 
   // PhiS
   PhiS = K->phiS;
+
+  // Theta
+  Double_t MRterm[2];
+  for(int h=0; h<2; h++)
+    MRterm[h] = TMath::Sqrt( vecH[h].M2() + vecR.Vect().Mag2() );
+  Theta = TMath::ACos(
+      ( MRterm[0] - MRterm[1] - Mh*Zeta ) / 
+      ( 2 * vecR.Vect().Mag() )  
+      );
 }
 
 ////////////////////////////////////////////////////////////////////

src/Dihadrons.h

@@ -11,7 +11,7 @@ class Dihadron {
     Dihadron() {};
     Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K);
     ~Dihadron() {};
-    Double_t Mh, MX, Z, PhPerp, Theta, PhiH, PhiR, PhiS;
+    Double_t Mh, MX, Z, PhPerp, Theta, PhiH, PhiR, PhiS, Zeta;
     Double_t hadZ[2];
   private:
   ClassDef(Dihadron,1);