avoid nan or inf when E and p^z are too close for particle samples

src/FSSW.cpp

@@ -1954,12 +1954,30 @@ int FSSW::sample_momemtum_from_a_fluid_cell(
 
             // assigned to the return variables
             pT = sqrt(pLab[1]*pLab[1] + pLab[2]*pLab[2]);
-            if (pT < 5.) {
-                phi = atan2(pLab[2], pLab[1]);
-                double y = 0.5*log((pLab[0] + pLab[3])/(pLab[0] - pLab[3]));
-                y_minus_eta_s = y - surf->eta;
-                return(1);
+            phi = atan2(pLab[2], pLab[1]);
+            double y = 0.5*log((pLab[0] + pLab[3])/(pLab[0] - pLab[3]));
+            if (std::isnan(y) || std::isinf(y)) {
+                // pLab[0] and pLab[3] are too close
+                // recompute pLab[3] with E, pT, and mass
+                pLab[3] = sqrt(pLab[0]*pLab[0] - pT*pT - mass*mass);
+                y = 0.5*log((pLab[0] + pLab[3])/(pLab[0] - pLab[3]));
+                if (std::isnan(y) || std::isinf(y)) {
+                    std::cout << "[Error]: sampled y is " << y << "!"
+                              << std::endl;
+                    std::cout << "umu = "
+                              << umu[0] << ", " << umu[1] << ", " << umu[2]
+                              << ", " << umu[3] << std::endl;
+                    std::cout << "pLRF = "
+                              << pLRF[0] << ", " << pLRF[1] << ", "
+                              << pLRF[2] << ", " << pLRF[3] << std::endl;
+                    std::cout << "pLab = "
+                              << pLab[0] << ", " << pLab[1] << ", "
+                              << pLab[2] << ", " << pLab[3] << std::endl;
+                    exit(1);
+                }
             }
+            y_minus_eta_s = y - surf->eta;
+            return(1);
         }
         tries++;
     }