trento/src/collider.cxx

0001 // TRENTO: Reduced Thickness Event-by-event Nuclear Topology
0002 // Copyright 2015 Jonah E. Bernhard, J. Scott Moreland
0003 // TRENTO3D: Three-dimensional extension of TRENTO by Weiyao Ke
0004 // MIT License
0005
0006 #include "collider.h"
0007
0008 #include <cmath>
0009 #include <string>
0010 #include <vector>
0011
0012 #include <boost/program_options/variables_map.hpp>
0013
0014 #include "fwd_decl.h"
0015 #include "nucleus.h"
0016 #include <iostream>
0017
0018 namespace trento {
0019
0020 namespace {
0021
0022 // Helper functions for Collider ctor.
0023
0024 // Create one nucleus from the configuration.
0025 NucleusPtr create_nucleus(const VarMap& var_map, std::size_t index) {
0026   const auto& species = var_map["projectile"]
0027                         .as<std::vector<std::string>>().at(index);
0028   const auto& nucleon_dmin = var_map["nucleon-min-dist"].as<double>();
0029   const auto& nucleon_width = var_map["nucleon-width"].as<double>();
0030   return Nucleus::create(species, nucleon_width, nucleon_dmin);
0031 }
0032
0033 // Determine the maximum impact parameter.  If the configuration contains a
0034 // non-negative value for bmax, use it; otherwise, fall back to the minimum-bias
0035 // default.
0036 double determine_bmax(const VarMap& var_map,
0037     const Nucleus& A, const Nucleus& B, const NucleonProfile& profile) {
0038   auto bmax = var_map["b-max"].as<double>();
0039   if (bmax < 0.)
0040     bmax = A.radius() + B.radius() + profile.max_impact();
0041   return bmax;
0042 }
0043
0044 // Determine the asymmetry parameter (Collider::asymmetry_) for a pair of
0045 // nuclei.  It's just rA/(rA+rB), falling back to 1/2 if both radii are zero
0046 // (i.e. for proton-proton).
0047 double determine_asym(const Nucleus& A, const Nucleus& B) {
0048   double rA = A.radius();
0049   double rB = B.radius();
0050   double sum = rA + rB;
0051   if (sum < 0.1)
0052     return 0.5;
0053   else
0054     return rA/sum;
0055 }
0056
0057 }  // unnamed namespace
0058
0059 // Lots of members to initialize...
0060 // Several helper functions are defined above.
0061 Collider::Collider(const VarMap& var_map)
0062     : nucleusA_(create_nucleus(var_map, 0)),
0063       nucleusB_(create_nucleus(var_map, 1)),
0064       nucleon_profile_(var_map),
0065       nevents_(var_map["number-events"].as<int>()),
0066       ntrys_(0),
0067       bmin_(var_map["b-min"].as<double>()),
0068       bmax_(determine_bmax(var_map, *nucleusA_, *nucleusB_, nucleon_profile_)),
0069       npartmin_(var_map["npart-min"].as<int>()),
0070       npartmax_(var_map["npart-max"].as<int>()),
0071       stotmin_(var_map["s-min"].as<double>()),
0072       stotmax_(var_map["s-max"].as<double>()),
0073       asymmetry_(determine_asym(*nucleusA_, *nucleusB_)),
0074       event_(var_map),
0075       output_(var_map),
0076       with_ncoll_(var_map["ncoll"].as<bool>())
0077 {
0078   // Constructor body begins here.
0079   // Set random seed if requested.
0080   auto seed = var_map["random-seed"].as<int64_t>();
0081   if (seed > 0)
0082     random::engine.seed(static_cast<random::Engine::result_type>(seed));
0083 }
0084
0085 // See header for explanation.
0086 Collider::~Collider() = default;
0087
0088 void Collider::run_events() {
0089   // The main event loop.
0090   for (int n = 0; n < nevents_; ++n) {
0091     if (n%1000 == 0 && n!=0) std::cout<< "# "<< n << " events generated" << std::endl;
0092     // Sampling the impact parameter also implicitly prepares the nuclei for
0093     // event computation, i.e. by sampling nucleon positions and participants.
0094
0095     // WK: an extra do-while loop, sample events, until it meets the Npart or
0096     // Entropy cut provided from command lines
0097     bool fullfil_Npart_cut=false, fullfil_Entropy_cut=false;
0098     double b;
0099     do{
0100         b = sample_impact_param();
0101         // Pass the prepared nuclei to the Event.  It computes the entropy profile
0102         // (thickness grid) and other event observables.
0103         event_.compute(*nucleusA_, *nucleusB_, nucleon_profile_);
0104         fullfil_Npart_cut = (npartmin_ < event_.npart())
0105                                 && (event_.npart() <= npartmax_);
0106         fullfil_Entropy_cut = (stotmin_ < event_.multiplicity())
0107                                 && (event_.multiplicity() <= stotmax_);
0108     }while( (!fullfil_Npart_cut) || (!fullfil_Entropy_cut) );
0109     // Write event data.
0110     output_(n, b, event_);
0111     records this_event;
0112     this_event.i = n;
0113     this_event.b = b;
0114     this_event.npart = event_.npart();
0115     this_event.mult = event_.multiplicity();
0116     all_records_.push_back(this_event);
0117   }
0118   double cross_section = nevents_*M_PI*(bmax_*bmax_ - bmin_*bmin_)/ntrys_;
0119   double cross_section_err = cross_section/std::sqrt(1.*nevents_);
0120   //std::cout << "# cross-section = " << cross_section
0121   //            << " +/- " << cross_section_err <<" [fm^2]" << std::endl;
0122 }
0123
0124 double Collider::sample_impact_param() {
0125   // Sample impact parameters until at least one nucleon-nucleon pair
0126   // participates.  The bool 'collision' keeps track -- it is effectively a
0127   // logical OR over all possible participant pairs.
0128   double b;
0129   bool collision = false;
0130
0131   do {
0132     // Sample b from P(b)db = 2*pi*b.
0133     b = std::sqrt(bmin_ * bmin_ + (bmax_ * bmax_ - bmin_ * bmin_) * random::canonical<double>());
0134
0135     // Offset each nucleus depending on the asymmetry parameter (see header).
0136     nucleusA_->sample_nucleons(asymmetry_ * b);
0137     nucleusB_->sample_nucleons((asymmetry_ - 1.) * b);
0138
0139     // Check each nucleon-nucleon pair.
0140     for (auto&& A : *nucleusA_) {
0141       for (auto&& B : *nucleusB_) {
0142         bool AB_collide = nucleon_profile_.participate(A, B);
0143
0144         if (with_ncoll_) {
0145             // WK: only init Ncoll and Ncoll density at the first binary collision:
0146             if (AB_collide && (!collision) ) event_.clear_TAB();
0147             // WK: to calculate binary collision denstiy, each collision
0148             // contribute independently its Tpp. Therefore, if one pair collide,
0149             // it calls the event object to accumulate Tpp to the Ncoll density
0150             // Ncoll density = Sum Tpp
0151             if (AB_collide) event_.accumulate_TAB(A, B, nucleon_profile_);
0152         }
0153
0154         // update collision flag
0155         collision = AB_collide || collision;
0156       }
0157     }
0158     ntrys_ ++;
0159   } while (!collision);
0160
0161   return b;
0162 }
0163
0164 }  // namespace trento