// No real positive solutions for exp(y1) found in Method 1 or 2.
//Method 3:
k1e = p1.e();
k1z = p1.z();
k2e = P.e() - p1.e();
k2z = P.z() - p1.z();
//This fixes transverse momenta
k1perp = sqrt(k1e*k1e-k1z*k1z);
k2perp = sqrt(k2e*k2e-k2z*k2z);
if (k1perp + k2perp > P.perp()){
std::cout<<"Implement Method 3."<<std::endl;
// Implement this if both Method 1 and 2 both fail.
} else {
//Method 4: We chose to set k2e = p2e and
double k2e = p2.e();
double k2z = p2.z();
double k1e = P.e() - p2.e();
double k1z = P.z() - p2.z();
//This fixes transverse momenta
double k1perp = sqrt(k1e*k1e-k1z*k1z);
double k2perp = sqrt(k2e*k2e-k2z*k2z);
if (k1perp + k2perp > P.perp()){
//Implement this if Method 1, 2 and 3 fail.
} else {
throw std::logic_error("No real positive solutions for exp(y) found in Method 1 or 2, and triangle inequality for transverse momenta is violated in Method 3 and 4.");
}
}
}
}
auto k1 = CLHEP::HepLorentzVector(k1x,k1y,k1z,k1e);
auto k2 = CLHEP::HepLorentzVector(k2x,k2y,k2z,k2e);