sPhenix code displayed by LXR master/​JETSCAPE/​external_packages/​trento/​test/​test_nucleus.cxx	Source navigation Diff markup Identifier search General search
	Version: master Revert   Change

File indexing completed on 2025-08-05 08:19:08

 // TRENTO: Reduced Thickness Event-by-event Nuclear Topology
 // Copyright 2015 Jonah E. Bernhard, J. Scott Moreland
 // MIT License
 
 #include "../src/nucleus.h"
 
 #include <algorithm>
 #include <cmath>
 #include <iomanip>
 #include <iterator>
 #include <map>
 
 #include "catch.hpp"
 #include "util.h"
 
 #include "../src/hdf5_utils.h"
 #include "../src/random.h"
 
 using namespace trento;
 
 TEST_CASE( "proton" ) {
   auto nucleus = Nucleus::create("p");
 
   CHECK( dynamic_cast<Proton*>(nucleus.get()) != nullptr );
 
   // proton contains one nucleon
   CHECK( std::distance(nucleus->begin(), nucleus->end()) == 1 );
   CHECK( std::distance(nucleus->cbegin(), nucleus->cend()) == 1 );
 
   // and has zero radius
   CHECK( nucleus->radius() == 0. );
 
   // sample position with random offset
   double offset = random::canonical<>();
   nucleus->sample_nucleons(offset);
   auto&& proton = *(nucleus->begin());
 
   // check correct position
   CHECK( proton.x() == offset );
   CHECK( proton.y() == 0. );
   CHECK( proton.z() == 0. );
 
   // not a participant initially
   CHECK( !proton.is_participant() );
 }
 
 TEST_CASE( "deuteron" ) {
   auto nucleus = Nucleus::create("d");
 
   CHECK( dynamic_cast<Deuteron*>(nucleus.get()) != nullptr );
 
   CHECK( std::distance(nucleus->begin(), nucleus->end()) == 2 );
   CHECK( std::distance(nucleus->cbegin(), nucleus->cend()) == 2 );
 
   nucleus->sample_nucleons(0);
 
   CHECK( nucleus->cbegin()->x() == -std::next(nucleus->cbegin())->x() );
   CHECK( nucleus->cbegin()->y() == -std::next(nucleus->cbegin())->y() );
   CHECK( nucleus->cbegin()->z() == -std::next(nucleus->cbegin())->z() );
 
   CHECK( nucleus->radius() == Approx(-std::log(.01)/(2*0.457)) );
 }
 
 TEST_CASE( "lead nucleus" ) {
   auto nucleus = Nucleus::create("Pb");
 
   CHECK( dynamic_cast<WoodsSaxonNucleus*>(nucleus.get()) != nullptr );
 
   constexpr int A = 208;
   CHECK( std::distance(nucleus->begin(), nucleus->end()) == A );
   CHECK( std::distance(nucleus->cbegin(), nucleus->cend()) == A );
 
   CHECK( nucleus->radius() > 6. );
 
   double offset = nucleus->radius() * random::canonical<>();
   nucleus->sample_nucleons(offset);
 
   // average nucleon position
   double x = 0., y = 0., z = 0.;
   for (const auto& nucleon : *nucleus) {
     x += nucleon.x();
     y += nucleon.y();
     z += nucleon.z();
   }
   x /= A;
   y /= A;
   z /= A;
   auto tolerance = .7;
   CHECK( std::abs(x - offset) < tolerance );
   CHECK( std::abs(y) < tolerance );
   CHECK( std::abs(z) < tolerance );
 
   // no initial participants
   bool initial_participants = std::any_of(
       nucleus->cbegin(), nucleus->cend(),
       [](decltype(*nucleus->cbegin())& n) {
         return n.is_participant();
       });
   CHECK( !initial_participants );
 }
 
 TEST_CASE( "copper nucleus" ) {
   auto nucleus = Nucleus::create("Cu");
   auto def_nucleus = Nucleus::create("Cu2");
 
   CHECK( dynamic_cast<WoodsSaxonNucleus*>(nucleus.get()) != nullptr );
   CHECK( dynamic_cast<DeformedWoodsSaxonNucleus*>(def_nucleus.get()) != nullptr );
 
   constexpr int A = 63;
   CHECK( std::distance(nucleus->begin(), nucleus->end()) == A );
   CHECK( std::distance(def_nucleus->cbegin(), def_nucleus->cend()) == A );
 
   CHECK( nucleus->radius() > 4. );
   CHECK( def_nucleus->radius() > nucleus->radius() );
 }
 
 TEST_CASE( "gold nucleus" ) {
   auto nucleus = Nucleus::create("Au");
   auto def_nucleus = Nucleus::create("Au2");
 
   CHECK( dynamic_cast<WoodsSaxonNucleus*>(nucleus.get()) != nullptr );
   CHECK( dynamic_cast<DeformedWoodsSaxonNucleus*>(def_nucleus.get()) != nullptr );
 
   constexpr int A = 197;
   CHECK( std::distance(nucleus->begin(), nucleus->end()) == A );
   CHECK( std::distance(def_nucleus->cbegin(), def_nucleus->cend()) == A );
 
   CHECK( nucleus->radius() > 6. );
   CHECK( def_nucleus->radius() > nucleus->radius() );
 }
 
 TEST_CASE( "uranium nucleus" ) {
   auto nucleus = Nucleus::create("U");
 
   CHECK( dynamic_cast<DeformedWoodsSaxonNucleus*>(nucleus.get()) != nullptr );
 
   constexpr int A = 238;
   CHECK( std::distance(nucleus->begin(), nucleus->end()) == A );
   CHECK( std::distance(nucleus->cbegin(), nucleus->cend()) == A );
 
   CHECK( nucleus->radius() > 8. );
 }
 
 #ifdef TRENTO_HDF5
 
 TEST_CASE( "manual nucleus" ) {
   std::vector<float> positions = {
      0.,  0., 0.,
      1.,  0., 0.,
     -1.,  0., 0.,
      3.,  0., 0.,
      0.,  2., 0.,
      0., -2., 0.,
   };
 
   const auto A = positions.size() / 3;
 
   temporary_path temp{".hdf5"};
 
   {
     auto dataspace = hdf5::make_dataspace(std::array<hsize_t, 3>{1, A, 3});
     const auto& datatype = hdf5::type<decltype(positions)::value_type>();
 
     H5::H5File file{temp.path.string(), H5F_ACC_EXCL};
 
     auto dataset = file.createDataSet("test", datatype, dataspace);
     dataset.write(positions.data(), datatype);
   }
 
   auto nucleus = Nucleus::create(temp.path.string());
 
   CHECK( dynamic_cast<ManualNucleus*>(nucleus.get()) != nullptr );
 
   CHECK( std::distance(nucleus->begin(), nucleus->end()) == A );
   CHECK( std::distance(nucleus->cbegin(), nucleus->cend()) == A );
 
   CHECK( nucleus->radius() == Approx(3.) );
 
   nucleus->sample_nucleons(0.);
 
   auto nucleon = nucleus->cbegin();
 
   CHECK( nucleon->x() == Approx(0.) );
   CHECK( nucleon->y() == Approx(0.) );
   CHECK( nucleon->z() == Approx(0.) );
 
   std::advance(nucleon, 1);
   CHECK( nucleon->x() == Approx(-std::next(nucleon)->x()) );
   CHECK( nucleon->y() == Approx(-std::next(nucleon)->y()) );
   CHECK( nucleon->z() == Approx(-std::next(nucleon)->z()) );
 
   CHECK( nucleon->x() == Approx(std::next(nucleon, 2)->x()/3) );
   CHECK( nucleon->y() == Approx(std::next(nucleon, 2)->y()/3) );
   CHECK( nucleon->z() == Approx(std::next(nucleon, 2)->z()/3) );
 
   std::advance(nucleon, 3);
   CHECK( nucleon->x() == Approx(-std::next(nucleon)->x()) );
   CHECK( nucleon->y() == Approx(-std::next(nucleon)->y()) );
   CHECK( nucleon->z() == Approx(-std::next(nucleon)->z()) );
 
   CHECK_THROWS_AS( Nucleus::create("nonexistent.hdf"), std::invalid_argument );
 }
 
 #endif  // TRENTO_HDF5
 
 TEST_CASE( "woods-saxon sampling" ) {
   int A = 200;
   double R = 6., a = .5;
   NucleusPtr nucleus{new WoodsSaxonNucleus{static_cast<std::size_t>(A), R, a}};
 
   // Test Woods-Saxon sampling.
   // This is honestly not a great test; while it does prove that the code
   // basically works as intended, it does not rigorously show that the generated
   // numbers are actually Woods-Saxon distributed.  Personally I am much more
   // convinced by plotting a histogram and visually comparing it to the smooth
   // curve.  The script 'plot-woods-saxon.py' in the 'woods-saxon' subdirectory
   // does this.
 
   // sample a bunch of nuclei and bin all the nucleon positions
   auto nevents = 5000;
   auto nsamples = nevents * A;
   std::map<int, int> hist{};
   for (auto i = 0; i < nevents; ++i) {
     nucleus->sample_nucleons(0.);
     for (const auto& nucleon : *nucleus) {
       auto x = nucleon.x();
       auto y = nucleon.y();
       auto z = nucleon.z();
       auto r = std::sqrt(x*x + y*y + z*z);
       ++hist[static_cast<int>(r)];
     }
   }
 
   // integrate the W-S dist from rmin to rmax
   auto integrate_woods_saxon = [R, a](double rmin, double rmax) -> double {
     auto nbins = static_cast<int>((rmax - rmin)/.001);
     auto dr = (rmax - rmin)/nbins;
     double result = 0.;
     for (auto n = 0; n <= nbins; ++n) {
       auto r = rmin + n*dr;
       auto f = r*r/(1. + std::exp((r - R)/a));
       if (n == 0 || n == nbins)
         f /= 2;
       result += f;
     }
     return result * dr;
   };
 
   double ws_norm = integrate_woods_saxon(0, hist.size());
 
   // check all histogram bins agree with numerical integration
   auto all_bins_correct = true;
   std::ostringstream bin_output{};
   bin_output << std::fixed << std::boolalpha
              << "rmin rmax  prob      cprob     ratio     pass\n";
   for (const auto& bin : hist) {
     auto rmin = bin.first;
     auto rmax = rmin + 1;
     auto prob = static_cast<double>(bin.second) / nsamples;
     auto correct_prob = integrate_woods_saxon(rmin, rmax) / ws_norm;
     bool within_tol = prob == Approx(correct_prob).epsilon(.1).margin(1e-4);
     if (!within_tol)
       all_bins_correct = false;
     bin_output << std::setw(4) << rmin << ' '
                << std::setw(4) << rmax << "  "
                << prob << "  "
                << correct_prob << "  "
                << prob/correct_prob << "  "
                << within_tol << '\n';
   }
   INFO( bin_output.str() );
   CHECK( all_bins_correct );
 }
 
 TEST_CASE( "deformed woods-saxon sampling" ) {
   // This is tough to test because the sampled positions are randomly rotated
   // and the z-coordinate is discarded.  However, the authors have visually
   // checked the results by histogramming the samples (not done here, but easy
   // to reproduce).
 
   // As a not-very-stringent test, check that the mean ellipticity of a deformed
   // W-S nucleus is significantly larger than a symmetric nucleus.
 
   std::size_t A = 200;
   double R = 6., a = .5, beta2 = .3, beta4 = .1;
 
   auto nucleus_sym = WoodsSaxonNucleus{A, R, a};
   auto nucleus_def = DeformedWoodsSaxonNucleus{A, R, a, beta2, beta4};
 
   auto mean_ecc2 = [](Nucleus& nucleus) {
     double e2 = 0.;
     int nevents = 500;
 
     for (int n = 0; n < nevents; ++n) {
       nucleus.sample_nucleons(0.);
       double numer = 0.;
       double denom = 0.;
       for (const auto& nucleon : nucleus) {
         double x2 = std::pow(nucleon.x(), 2);
         double y2 = std::pow(nucleon.y(), 2);
         numer += x2 - y2;
         denom += x2 + y2;
       }
 
       e2 += std::fabs(numer) / denom;
     }
 
     return e2 / nevents;
   };
 
   CHECK( mean_ecc2(nucleus_def) > 2.*mean_ecc2(nucleus_sym) );
 }
 
 TEST_CASE( "nuclear radius" ) {
   constexpr auto R = 5., a = .5;
   constexpr auto radius = R + 3*a;
 
   WoodsSaxonNucleus nucleus{100, R, a};
   DeformedWoodsSaxonNucleus dummy_def_nucleus{100, R, a, 0., 0.};
 
   CHECK( nucleus.radius() == Approx(radius) );
   CHECK( dummy_def_nucleus.radius() == Approx(radius) );
 
   DeformedWoodsSaxonNucleus def_nucleus{100, R, a, .2, .1};
   CHECK( def_nucleus.radius() > radius );
 }
 
 TEST_CASE( "minimum distance" ) {
   for (const auto& species : {"Pb", "U"}) {
     for (auto repeat = 0; repeat < 10; ++repeat) {
       const auto target_dmin = .2 + .4*random::canonical<>();
       auto nucleus = Nucleus::create("Pb", target_dmin);
       nucleus->sample_nucleons(10*random::canonical<>());
       auto dminsq = 100.;
       for (auto n1 = nucleus->cbegin(); n1 != nucleus->cend(); ++n1) {
         for (auto n2 = n1 + 1; n2 != nucleus->cend(); ++n2) {
           auto dx = n1->x() - n2->x();
           auto dy = n1->y() - n2->y();
           auto dz = n1->z() - n2->z();
           dminsq = std::fmin(dx*dx + dy*dy + dz*dz, dminsq);
         }
       }
       auto dmin = std::sqrt(dminsq);
       INFO( "species " << species << " repeat " << repeat );
       CHECK( dmin >= target_dmin );
     }
   }
 }
 
 TEST_CASE( "unknown nucleus species" ) {
   CHECK_THROWS_AS( Nucleus::create("hello"), std::invalid_argument );
 }

This page was automatically generated by the 2.3.7 LXR engine. The LXR team