trento/test/test_nucleon.cxx

0001 // TRENTO: Reduced Thickness Event-by-event Nuclear Topology
0002 // Copyright 2015 Jonah E. Bernhard, J. Scott Moreland
0003 // MIT License
0004
0005 #include "../src/nucleon.h"
0006
0007 #include <cmath>
0008
0009 #include "catch.hpp"
0010 #include "util.h"
0011
0012 #include "../src/nucleus.h"
0013 #include "../src/random.h"
0014
0015 using namespace trento;
0016
0017 TEST_CASE( "nucleon" ) {
0018   auto fluct = 1. + .5*random::canonical<>();
0019   auto xsec = 4. + 3.*random::canonical<>();
0020   auto width = .5 + .2*random::canonical<>();
0021   auto wsq = width*width;
0022
0023   auto var_map = make_var_map({
0024       {"fluctuation",   fluct},
0025       {"cross-section", xsec},
0026       {"nucleon-width", width},
0027   });
0028
0029   NucleonProfile profile{var_map};
0030   profile.fluctuate();
0031
0032   // truncation radius
0033   auto R = profile.radius();
0034   CHECK( R == Approx(5*width) );
0035
0036   // thickness function
0037   // check relative to zero
0038   auto tzero = profile.thickness(0.);
0039   CHECK( profile.thickness(wsq) == Approx(tzero*std::exp(-.5)) );
0040
0041   // random point inside radius
0042   auto dsq = std::pow(R*random::canonical<>(), 2);
0043   CHECK( profile.thickness(dsq) == Approx(tzero*std::exp(-.5*dsq/wsq)).epsilon(1e-5).margin(1e-5) );
0044
0045   // random point outside radius
0046   dsq = std::pow(R*(1+random::canonical<>()), 2);
0047   CHECK( profile.thickness(dsq) == 0. );
0048
0049   // fluctuations
0050   // just check they have unit mean -- the rest is handled by the C++ impl.
0051   auto total = 0.;
0052   auto n = 1e6;
0053   for (auto i = 0; i < static_cast<int>(n); ++i) {
0054     profile.fluctuate();
0055     total += profile.thickness(0.) * (2*M_PI*wsq);
0056   }
0057
0058   auto mean = total/n;
0059   CHECK( mean == Approx(1.).epsilon(.003) );
0060
0061   // must use a Nucleus to set Nucleon position
0062   // Proton conveniently sets a deterministic position
0063   // a mock class would be better but this works fine
0064   Proton A{}, B{};
0065   A.sample_nucleons(0.);
0066   B.sample_nucleons(0.);
0067   auto& nA = *A.begin();
0068   auto& nB = *B.begin();
0069   CHECK( nA.x() == 0. );
0070   CHECK( nA.y() == 0. );
0071   CHECK( nA.z() == 0. );
0072   CHECK( !nA.is_participant() );
0073
0074   // wait until the nucleons participate
0075   while (!profile.participate(nA, nB)) {}
0076   CHECK( nA.is_participant() );
0077   CHECK( nB.is_participant() );
0078
0079   // resampling nucleons resets participant state
0080   A.sample_nucleons(0.);
0081   CHECK( !nA.is_participant() );
0082
0083   // test cross section
0084   // min-bias impact params
0085   auto bmax = profile.max_impact();
0086   CHECK( bmax == Approx(6*width) );
0087
0088   auto nev = 1e6;
0089   auto count = 0;
0090   for (auto i = 0; i < static_cast<int>(nev); ++i) {
0091     auto b = bmax * std::sqrt(random::canonical<>());
0092     A.sample_nucleons(.5*b);
0093     B.sample_nucleons(-.5*b);
0094     if (profile.participate(nA, nB))
0095       ++count;
0096   }
0097
0098   auto xsec_mc = M_PI*bmax*bmax * static_cast<double>(count)/nev;
0099
0100   // precision is better than this, but let's be conservative
0101   CHECK( xsec_mc == Approx(xsec).epsilon(.02) );
0102
0103   // impact larger than max should never participate
0104   auto b = bmax + random::canonical<>();
0105   A.sample_nucleons(.5*b);
0106   B.sample_nucleons(-.5*b);
0107   CHECK( !profile.participate(nA, nB) );
0108
0109   // very large fluctuation parameters mean no fluctuations
0110   auto no_fluct_var_map = make_var_map({
0111       {"fluctuation",   1e12},
0112       {"cross-section", xsec},
0113       {"nucleon-width", width},
0114   });
0115
0116   NucleonProfile no_fluct_profile{no_fluct_var_map};
0117   no_fluct_profile.fluctuate();
0118   CHECK( no_fluct_profile.thickness(0) == Approx(1/(2*M_PI*wsq)) );
0119
0120   CHECK_THROWS_AS([]() {
0121     // nucleon width too small
0122     auto bad_var_map = make_var_map({
0123         {"fluctuation",   1.},
0124         {"cross-section", 5.},
0125         {"nucleon-width", .1},
0126     });
0127     NucleonProfile bad_profile{bad_var_map};
0128   }(),
0129   std::domain_error);
0130 }