src/framework/InitialState.cc

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015
0016 #include "InitialState.h"
0017 #include "JetScapeWriter.h"
0018 #include <iostream>
0019
0020 namespace Jetscape {
0021
0022 InitialState::~InitialState() {}
0023
0024 void InitialState::Init() {
0025   JetScapeModuleBase::Init();
0026
0027   JSINFO << "Initialize InitialState ... " << GetId() << " ...";
0028
0029   grid_max_x_ = GetXMLElementDouble({"IS", "grid_max_x"});
0030   grid_max_y_ = GetXMLElementDouble({"IS", "grid_max_y"});
0031   grid_max_z_ = GetXMLElementDouble({"IS", "grid_max_z"});
0032   grid_step_x_ = GetXMLElementDouble({"IS", "grid_step_x"});
0033   grid_step_y_ = GetXMLElementDouble({"IS", "grid_step_y"});
0034   grid_step_z_ = GetXMLElementDouble({"IS", "grid_step_z"});
0035   JSINFO << "x range for bulk evolution = [" << -grid_max_x_ << ", "
0036          << grid_max_x_ << "]";
0037
0038   InitTask();
0039
0040   JetScapeTask::InitTasks();
0041 }
0042
0043 void InitialState::Exec() {
0044   // Do whatever is needed to figure out the internal temp...
0045 }
0046
0047 void InitialState::Clear() {}
0048
0049 void InitialState::Write(weak_ptr<JetScapeWriter> w) {
0050   //Write out the original vertex so the writer can keep track of it...
0051   // auto f = w.lock();
0052   // if ( f ) f->Write(make_shared<Vertex>(initialVtx));
0053 }
0054
0055 void InitialState::CollectHeader(weak_ptr<JetScapeWriter> w) {
0056   auto f = w.lock();
0057   if (f) {
0058     auto &header = f->GetHeader();
0059     header.SetNpart(GetNpart());
0060     header.SetNcoll(GetNcoll());
0061     header.SetTotalEntropy(GetTotalEntropy());
0062   }
0063 }
0064
0065 std::tuple<double, double, double> InitialState::CoordFromIdx(int idx) {
0066   int nx = GetXSize();
0067   int ny = GetYSize();
0068   int nz = GetZSize();
0069
0070   int page = idx / (nx * ny);
0071   int row = (idx - page * nx * ny) / nx;
0072   int col = idx - page * nx * ny - row * nx;
0073
0074   return std::make_tuple(-grid_max_x_ + col * grid_step_x_,
0075                          -grid_max_y_ + row * grid_step_y_,
0076                          -grid_max_z_ + page * grid_step_z_);
0077 }
0078
0079
0080 void InitialState::SampleABinaryCollisionPoint(double &x, double &y) {
0081   if (num_of_binary_collisions_.size() == 0) {
0082     JSWARN << "num_of_binary_collisions is empty, setting the starting "
0083               "location to 0. Make sure to add e.g. trento before PythiaGun.";
0084   } else {
0085     std::discrete_distribution<> dist(
0086         begin(num_of_binary_collisions_),
0087         end(num_of_binary_collisions_)); // Create the distribution
0088     // Now generate values
0089     auto idx = dist(*GetMt19937Generator());
0090     auto coord = CoordFromIdx(idx);
0091     x = std::get<0>(coord);
0092     y = std::get<1>(coord);
0093   }
0094 }
0095
0096 } // end namespace Jetscape