Merge pull request #74 from JETSCAPE/Passing-QCDStringList-to-MUSIC

Passing qcd string list to music

Author: Hendrik1704

.github/workflows/test-events-Pythia-Isr-Hadron.yaml

@@ -55,7 +55,7 @@ jobs:
       uses: actions/checkout@v4
       with:
         repository: JETSCAPE/TEST-EXAMPLES
-        ref: main
+        ref: Passing-QCDStringList-to-MUSIC
         path: TEST-EXAMPLES
 
     - name: Run Pythia ISR Hadron File Tests

config/jetscape_main.xml

@@ -400,7 +400,9 @@
       <!--higher order transport coefficients-->
       <Include_second_order_terms>1</Include_second_order_terms>
       <!--freeze-out parameters-->
+      <use_eps_for_freeze_out>0</use_eps_for_freeze_out>
       <freezeout_temperature>0.15</freezeout_temperature>
+      <eps_switch>0.3</eps_switch>
       <!--perform CooperFrye Freezout in MUSIC, not needed when iSS module is used-->
       <Perform_CooperFrye_Freezeout>0</Perform_CooperFrye_Freezeout>
       <!--source term parameters for InitialProfile 13-->

external_packages/get_3dglauber.sh

@@ -14,7 +14,7 @@
 ##############################################################################
 
 folderName="3dMCGlauber"
-commitHash="056a6d03bf6fb097b14d82874cebf6cce2fa0e5a" # for xscape 1.1, 1.1.1
+commitHash="b54d997a8b35" # for xscape 1.2
 # download the code package
 rm -fr $folderName
 git clone https://github.com/chunshen1987/3dMCGlauber.git --branch JETSCAPE $folderName

src/framework/HardProcess.cc

@@ -123,9 +123,12 @@ void HardProcess::CollectHeader(weak_ptr<JetScapeWriter> w) {
     header.SetSigmaErr(GetSigmaErr());
     header.SetPtHat(GetPtHat());
     header.SetEventWeight(GetEventWeight());
-    header.SetVertexX(hp_list[0]->x_in().x());
-    header.SetVertexY(hp_list[0]->x_in().y());
-    header.SetVertexZ(hp_list[0]->x_in().z());
+
+    if (hp_list.size() > 0) {
+      header.SetVertexX(hp_list[0]->x_in().x());
+      header.SetVertexY(hp_list[0]->x_in().y());
+      header.SetVertexZ(hp_list[0]->x_in().z());
+    }
   }
 }
 

src/framework/InitialState.h

@@ -131,17 +131,17 @@ class InitialState : public JetScapeModuleBase {
   std::tuple<double, double, double> CoordFromIdx(int idx);
   virtual void SampleABinaryCollisionPoint(double &t, double &x,
                                            double &y, double &z);
-                                           
+
   virtual double Get_total_nucleon_density_lab(double t, double x,
                                                double y, double z) {
       return 0.0;
   }
-                                           
+
   virtual double Get_target_nucleon_density_lab(double t, double x,
                                                 double y, double z) {
       return 0.0;
   }
-                                                
+
   virtual double Get_projectile_nucleon_density_lab(double t, double x,
                                                     double y, double z) {
       return 0.0;
@@ -153,7 +153,6 @@ class InitialState : public JetScapeModuleBase {
   virtual void OutputHardPartonMomentum(double E, double px, double py, double pz,
                                         int direction, double P_A);
   virtual void ClearHardPartonMomentum();
-                                   
 
   virtual std::vector<double> Get_projectile_nucleon_z_lab();
   virtual std::vector<double> Get_target_nucleon_z_lab();
@@ -169,6 +168,9 @@ class InitialState : public JetScapeModuleBase {
 
   virtual void GetHardPartonPosAndMomentumProj();
   virtual void GetHardPartonPosAndMomentumTarg();
+  virtual std::vector< std::vector<double> > GetQCDStringList() {
+    return(QCDStringList_);
+}
 
   virtual void GenerateStrings();
   /**  @return The maximum value of coordinate "x" in the nuclear profile of a nucleus.
@@ -237,6 +239,7 @@ class InitialState : public JetScapeModuleBase {
       @sa Function CoordFromIdx(int idx) for the mapping.
    */
   std::vector<double> num_of_binary_collisions_;
+  std::vector< std::vector<double> > QCDStringList_;
   // the above should be private. Only Adding getters for now to not break other people's code
 
   /**  @return The initial state entropy density distribution.

src/hydro/MusicWrapper.cc

@@ -111,9 +111,9 @@ void MpiMusic::InitializeHydro(Parameter parameter_list) {
     exit(1);
   }
 
-  int InitialProfile = (
+  initialProfile_ = (
       GetXMLElementInt({"Hydro", "MUSIC", "InitialProfile"}));
-  music_hydro_ptr->set_parameter("Initial_profile", InitialProfile);
+  music_hydro_ptr->set_parameter("Initial_profile", initialProfile_);
   double string_source_sigma_x = (
         GetXMLElementDouble({"Hydro", "MUSIC", "string_source_sigma_x"}));
   music_hydro_ptr->set_parameter("string_source_sigma_x",
@@ -207,6 +207,10 @@ void MpiMusic::InitializeHydro(Parameter parameter_list) {
     music_hydro_ptr->set_parameter("Include_second_order_terms", 1);
   }
 
+  int use_eps_for_freeze_out = GetXMLElementInt(
+    {"Hydro", "MUSIC", "use_eps_for_freeze_out"});
+  music_hydro_ptr->set_parameter("use_eps_for_freeze_out", 
+    use_eps_for_freeze_out);
   freezeout_temperature =
       GetXMLElementDouble({"Hydro", "MUSIC", "freezeout_temperature"});
   if (freezeout_temperature > 0.05) {
@@ -216,7 +220,10 @@ void MpiMusic::InitializeHydro(Parameter parameter_list) {
            << freezeout_temperature << " GeV!";
     exit(1);
   }
-
+  double eps_switch =
+      GetXMLElementDouble({"Hydro", "MUSIC", "eps_switch"});
+  music_hydro_ptr->set_parameter("eps_switch", eps_switch);
+  
   music_hydro_ptr->check_parameters();
   music_hydro_ptr->add_hydro_source_terms(hydro_source_terms_ptr);
 }
@@ -231,22 +238,28 @@ void MpiMusic::InitializeHydroEnergyProfile() {
   int nz = ini->GetZSize();
 
   // need further improvement to accept multiple source term objects
-  music_hydro_ptr->generate_hydro_source_terms();
-
-  if (pre_eq_ptr == nullptr) {
-    JSWARN << "Missing the pre-equilibrium module ...";
-    exit(1);
-  }
-  double tau0 = pre_eq_ptr->GetPreequilibriumEndTime();
-  JSINFO << "hydro initial time  tau0 = " << tau0 << " fm";
-  music_hydro_ptr->initialize_hydro_from_jetscape_preequilibrium_vectors(
+  if (initialProfile_ == 13 || initialProfile_ == 131) {
+    music_hydro_ptr->generate_hydro_source_terms(ini->GetQCDStringList());
+    music_hydro_ptr->initialize_hydro_xscape();
+  } else {
+    music_hydro_ptr->generate_hydro_source_terms();
+    double tau0 = pre_eq_ptr->GetPreequilibriumEndTime();
+    JSINFO << "hydro initial time  tau0 = " << tau0 << " fm";
+    music_hydro_ptr->initialize_hydro_from_jetscape_preequilibrium_vectors(
       tau0,
       dx, dz, z_max, nz, pre_eq_ptr->e_, pre_eq_ptr->P_,
       pre_eq_ptr->utau_, pre_eq_ptr->ux_, pre_eq_ptr->uy_, pre_eq_ptr->ueta_,
       pre_eq_ptr->pi00_, pre_eq_ptr->pi01_, pre_eq_ptr->pi02_,
       pre_eq_ptr->pi03_, pre_eq_ptr->pi11_, pre_eq_ptr->pi12_,
       pre_eq_ptr->pi13_, pre_eq_ptr->pi22_, pre_eq_ptr->pi23_,
       pre_eq_ptr->pi33_, pre_eq_ptr->bulk_Pi_);
+  }
+
+  if (pre_eq_ptr == nullptr) {
+    JSWARN << "Missing the pre-equilibrium module ...";
+    exit(1);
+  }
+
   JSINFO << "initial density profile dx = " << dx << " fm";
   hydro_status = INITIALIZED;
   JSINFO << "number of source terms: "

src/hydro/MusicWrapper.h

@@ -91,6 +91,8 @@ class MpiMusic : public FluidDynamics {
   bool has_source_terms;
   std::shared_ptr<HydroSourceJETSCAPE> hydro_source_terms_ptr;
 
+  int initialProfile_;
+
   // Allows the registration of the module so that it is available to be
   // used by the Jetscape framework.
   static RegisterJetScapeModule<MpiMusic> reg;

src/initialstate/MCGlauberWrapper.cc

@@ -34,13 +34,12 @@ MCGlauberWrapper::MCGlauberWrapper() {
 
 
 void MCGlauberWrapper::InitTask() {
-    parameter_list_.read_in_parameters_from_file("mcglauber.input");
-    //int ran_seed = parameter_list_.get_seed();
     auto ran_seed = (*GetMt19937Generator())();
-    auto gamma_beta = parameter_list_.get_tau_form_fluct_gamma_beta();
 
+    int argc = 0;
+    char* argv[1];
     mc_gen_ = std::shared_ptr<MCGlb::EventGenerator>(
-              new MCGlb::EventGenerator("mcglauber.input", ran_seed));
+        new MCGlb::EventGenerator("mcglauber.input", argc, argv, ran_seed));
 
     // overwrite input options
     double para_temp;
@@ -65,21 +64,19 @@ void MCGlauberWrapper::InitTask() {
     para_temp = (
         GetXMLElementDouble({"IS", "MCGlauber", "ylossParam4var"}));
     mc_gen_->set_parameter("ylossParam4var", para_temp);
-    para_temp = (
-        GetXMLElementDouble({"IS", "MCGlauber", "remnant_energy_loss_fraction"}));
+    para_temp = (GetXMLElementDouble({"IS", "MCGlauber",
+                                      "remnant_energy_loss_fraction"}));
     mc_gen_->set_parameter("remnant_energy_loss_fraction", para_temp);
-
-    // re-generate mc pointer
-    mc_gen_->New_Para_pointer(ran_seed);
 }
 
+
 void MCGlauberWrapper::ClearTask() {
     VERBOSE(1) << "clear initial condition vectors";
     binary_collision_t_.clear();
     binary_collision_x_.clear();
     binary_collision_y_.clear();
     binary_collision_z_.clear();
-
+    QCDStringList_.clear();
 }
 
 
@@ -115,7 +112,7 @@ void MCGlauberWrapper::ExecuteTask() {
 
 void MCGlauberWrapper::SampleABinaryCollisionPoint(
         double &t, double &x, double &y, double &z) {
-    int rand_idx = (*rand_int_ptr_)(*GetMt19937Generator());
+    const int rand_idx = (*rand_int_ptr_)(*GetMt19937Generator());
     t = binary_collision_t_[rand_idx];
     x = binary_collision_x_[rand_idx];
     y = binary_collision_y_[rand_idx];
@@ -206,54 +203,74 @@ void MCGlauberWrapper::OutputHardPartonMomentum(double E, double px, double py,
                              << " targ_parton_e_ " << targ_parton_e_ 
                              << " targ_parton_pz_ " << targ_parton_pz_;
 
-    if(targ_parton_e_ >= 0.95 * P_A || proj_parton_e_ >= 0.95 * P_A ){
-        throw std::runtime_error("Energy to subtract from 3DMCGlauber >= 0.95 * P_A "+std::to_string(0.95 * P_A)+". Turn on Verbose for more info.");
+    if (targ_parton_e_ >= 0.95 * P_A || proj_parton_e_ >= 0.95 * P_A ) {
+        throw std::runtime_error(
+            "Energy to subtract from 3DMCGlauber >= 0.95 * P_A "
+            + std::to_string(0.95 * P_A)+". Turn on Verbose for more info.");
     }
 
 }
 
+
 std::vector<double> MCGlauberWrapper::Get_quarks_pos_proj_lab() {
     // get the x, y, z of the three valence quarks of colliding projectile
     // The fourth parton is the soft ball
-    // 3DGlauber attributes the remaining energy and momentum carried by the sea quarks and gluons to a soft gluon cloud
+    // 3DGlauber attributes the remaining energy and momentum carried by the
+    // sea quarks and gluons to a soft gluon cloud
     // Output formulation is (x,y,z, x,y,z, x,y,z, x,y,z)
-    mc_gen_->GetHardPos(hard_parton_t_, hard_parton_x_, hard_parton_y_, 
+    mc_gen_->GetHardPos(hard_parton_t_, hard_parton_x_, hard_parton_y_,
                         hard_parton_z_);
     return(mc_gen_->GetQuarkPosProj());
 }
 
+
 std::vector<double> MCGlauberWrapper::Get_quarks_pos_targ_lab() {
     // get the x, y, z of the three valence quarks of colliding target
     // The fourth parton is the soft ball
-    // 3DGlauber attributes the remaining energy and momentum carried by the sea quarks and gluons to a soft gluon cloud
+    // 3DGlauber attributes the remaining energy and momentum carried by the
+    // sea quarks and gluons to a soft gluon cloud
     // Output formulation is (x,y,z, x,y,z, x,y,z, x,y,z)
-    mc_gen_->GetHardPos(hard_parton_t_, hard_parton_x_, hard_parton_y_, 
+    mc_gen_->GetHardPos(hard_parton_t_, hard_parton_x_, hard_parton_y_,
                         hard_parton_z_);
     return(mc_gen_->GetQuarkPosTarg());
 }
 
+
 std::vector<double> MCGlauberWrapper::Get_remnant_proj() {
     // get the fout-momentum (E, px, py, pz) of the remnant in projectile
     return(mc_gen_->GetRemMom_Proj());
 }
 
+
 std::vector<double> MCGlauberWrapper::Get_remnant_targ() {
-    // get the fout-momentum (E, px, py, pz) of the remnant in target 
+    // get the fout-momentum (E, px, py, pz) of the remnant in target
     return(mc_gen_->GetRemMom_Targ());
 }
 
+
 void MCGlauberWrapper::GetHardPartonPosAndMomentumProj() {
-    mc_gen_->GetMomandPos_Proj(hard_parton_t_, hard_parton_x_, hard_parton_y_, 
-                               hard_parton_z_, proj_parton_e_, proj_parton_px_, 
+    mc_gen_->GetMomandPos_Proj(hard_parton_t_, hard_parton_x_, hard_parton_y_,
+                               hard_parton_z_, proj_parton_e_, proj_parton_px_,
                                proj_parton_py_, proj_parton_pz_);
 }
 
+
 void MCGlauberWrapper::GetHardPartonPosAndMomentumTarg() {
-    mc_gen_->GetMomandPos_Targ(hard_parton_t_, hard_parton_x_, hard_parton_y_, 
-                               hard_parton_z_, targ_parton_e_, targ_parton_px_, 
+    mc_gen_->GetMomandPos_Targ(hard_parton_t_, hard_parton_x_, hard_parton_y_,
+                               hard_parton_z_, targ_parton_e_, targ_parton_px_,
                                targ_parton_py_, targ_parton_pz_);
 }
 
 void MCGlauberWrapper::GenerateStrings() {
-    mc_gen_->generate_strings(); // generate 3D Glauber for MUSIC
+    // generate strings from 3D Glauber for MUSIC
+    QCDStringList_.clear();
+    auto stringList = mc_gen_->generate_strings();
+    for (auto string_i: stringList) {
+        std::vector<double> string_temp;
+        for (auto ii: string_i) {
+            string_temp.push_back(static_cast<double>(ii));
+        }
+        QCDStringList_.push_back(string_temp);
+    }
+    JSINFO << "Produced " << QCDStringList_.size() << " strings.";
 }

src/initialstate/MCGlauberWrapper.h

@@ -35,7 +35,7 @@ class MCGlauberWrapper : public Jetscape::InitialState {
 public:
   MCGlauberWrapper();
   ~MCGlauberWrapper() {}
-  
+
   /** Reads the input parameters from the XML file under the tag  <IS>. Calls InitTask(); This explicit call of InitTask() can be used for actual initialization of modules such as @a Trento if attached as a @a polymorphic class. It also initializes the tasks within the current module.
       @sa Read about @a polymorphism in C++.
    */
@@ -71,7 +71,7 @@ class MCGlauberWrapper : public Jetscape::InitialState {
   void OutputHardCollisionPosition(double t, double x, double y, 
                                    double z);
   void OutputHardPartonMomentum(double E, double px, double py, double pz,
-                                                           int direction, double P_A);
+                                int direction, double P_A);
   void ClearHardPartonMomentum();
   void GetHardPartonPosAndMomentumProj();
   void GetHardPartonPosAndMomentumTarg();
@@ -82,19 +82,24 @@ class MCGlauberWrapper : public Jetscape::InitialState {
   std::vector<double> Get_remnant_proj();
   std::vector<double> Get_remnant_targ();
   void GenerateStrings();
-  
-  std::shared_ptr<InitialState> ini; 
-  
+  std::vector< std::vector<double> > GetQCDStringList() {
+    return(QCDStringList_);
+}
+
+std::shared_ptr<InitialState> ini;
+
 private:
   std::shared_ptr<MCGlb::EventGenerator> mc_gen_;
   std::vector<double> binary_collision_t_;
   std::vector<double> binary_collision_x_;
   std::vector<double> binary_collision_y_;
   std::vector<double> binary_collision_z_;
+  std::vector< std::vector<double> > QCDStringList_;
   double hard_parton_x_, hard_parton_y_, hard_parton_z_, hard_parton_t_;
-  double targ_parton_px_ = 0.0, targ_parton_py_ = 0.0, targ_parton_pz_ = 0.0, targ_parton_e_ = 0.0;
-  double proj_parton_px_ = 0.0, proj_parton_py_ = 0.0, proj_parton_pz_ = 0.0, proj_parton_e_ = 0.0;
-  MCGlb::Parameters parameter_list_;
+  double targ_parton_px_ = 0.0, targ_parton_py_ = 0.0;
+  double targ_parton_pz_ = 0.0, targ_parton_e_ = 0.0;
+  double proj_parton_px_ = 0.0, proj_parton_py_ = 0.0;
+  double proj_parton_pz_ = 0.0, proj_parton_e_ = 0.0;
   std::shared_ptr<MCGlb::RandomUtil::Random> ran_gen_ptr_;
   std::shared_ptr<std::uniform_int_distribution<int>> rand_int_ptr_;
   int ncoll_ = -1;