Bug fixes, C++17, softMomentumCutoff #168

Bug fixes, C++17, softMomentumCutoff #168
From JETSCAPE 3.6

CMakeLists.txt

@@ -102,7 +102,7 @@ endif (USE_SMASH)
 ### Compiler & Linker Flags ###
 ###############################
 message("Compiler and Linker flags ...")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC -pipe -Wall -std=c++11")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC -pipe -Wall -std=c++17")
 ## can turn off some warnings
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-reorder -Wno-unused-variable -Wno-inconsistent-missing-override -Wno-unused-private-field")
 ## Then set the build type specific options. These will be automatically appended.

config/jetscape_main.xml

@@ -116,6 +116,8 @@
       <pTHatMin>100</pTHatMin>
       <pTHatMax>120</pTHatMax>
       <eCM>5020</eCM>
+      <!-- The soft momentum cutoff in GeV removes partons from the pythia gun. The default matches the PP19 setup. -->
+      <softMomentumCutoff>2.0</softMomentumCutoff>
       <useHybridHad>0</useHybridHad>
       <!-- You can add any number of additional lines to initialize pythia here -->
       <!-- Note that if the tag exists it cannot be empty (tinyxml produces a segfault) -->

config/jetscape_user_FSFileTest.xml

@@ -26,7 +26,6 @@
       <PGun>
         <name>PGun</name>
         <pT>100</pT>
-        <useHybridHad>0</useHybridHad>
       </PGun>
   </Hard>
 

config/jetscape_user_HydroFile.xml

@@ -26,7 +26,6 @@
       <PGun>
         <name>PGun</name>
         <pT>100</pT>
-        <useHybridHad>0</useHybridHad>
       </PGun>
   </Hard>
 

config/jetscape_user_pbpb_grid.xml

@@ -1,11 +1,11 @@
 <?xml version="1.0"?>
 
 <jetscape>
-  
+
   <nEvents> 5 </nEvents>
   <setReuseHydro> true </setReuseHydro>
   <nReuseHydro> 10 </nReuseHydro>
-  
+
   <JetScapeWriterAscii> on </JetScapeWriterAscii>
   <outputFilename>test_out</outputFilename>
   <vlevel> 0 </vlevel>
@@ -19,7 +19,7 @@
     <!-- initial_profile_path  /wsu/home/groups/maj-shen/JETSCAPEDataFile/HydroProfiles/5020GeV_0-10  initial_profile_path -->
      <initial_profile_path>../examples/test_hydro_files</initial_profile_path>
   </IS>
-  
+
   <!-- Hard Process -->
   <Hard>
     <PythiaGun>
@@ -30,24 +30,23 @@
       <!-- You can add any number of additional lines to initialize pythia here -->
       <!-- Note that if the tag exists it cannot be empty (tinyxml produces a segfault) -->
       <!-- EPS09 Avg Nuclear PDF for nucleon in Pb nucleus -->
-      <LinesToRead> 
+      <LinesToRead>
 	HardQCD:all = on
         PDF:useHardNPDFA=on
 	PDF:useHardNPDFB=on
-	PDF:nPDFSetA=1         
+	PDF:nPDFSetA=1
 	PDF:nPDFSetB=1
 	PDF:nPDFBeamA = 100822080
 	PDF:nPDFBeamB = 100822080
       </LinesToRead>
-      <useHybridHad>0</useHybridHad>
     </PythiaGun>
   </Hard>
-  
+
   <!--Preequilibrium Dynamics Module -->
   <Preequilibrium>
     <NullPreDynamics> </NullPreDynamics>
   </Preequilibrium>
-  
+
   <!-- Hydro  Module  -->
   <Hydro>
     <hydro_from_file>
@@ -57,7 +56,7 @@
       <hydro_files_folder>../examples/test_hydro_files</hydro_files_folder>
     </hydro_from_file>
   </Hydro>
-  
+
   <!--Eloss Modules -->
   <Eloss>
 
@@ -108,7 +107,7 @@
       <run_alphas>1</run_alphas>   <!-- 0 for fixed alpha_s and 1 for running alpha_s -->
     </Lbt>
   </Eloss>
-  
+
   <!-- Jet Hadronization Module -->
   <JetHadronization>
     <name>colorless</name>
@@ -118,5 +117,5 @@
       111:mayDecay = off
     </LinesToRead>
   </JetHadronization>
-  
+
 </jetscape>

config/jetscape_user_twostagehydro.xml

@@ -6,7 +6,7 @@
   <setReuseHydro> true </setReuseHydro>
   <nReuseHydro> 5 </nReuseHydro>
   <vlevel> 0 </vlevel>
-  
+
   <!-- Jetscape Writer -->
   <JetScapeWriterAscii> on </JetScapeWriterAscii>
 
@@ -20,7 +20,6 @@
       <PGun>
         <name>PGun</name>
         <pT>100</pT>
-        <useHybridHad>0</useHybridHad>
       </PGun>
   </Hard>
 
@@ -35,7 +34,7 @@
       <name>MUSIC_1</name>
     </MUSIC>
   </Hydro>
-  
+
   <!-- Create liquifier -->
   <Liquefier>
       <!-- CausalLiquefier -->
@@ -51,7 +50,7 @@
           <width_delta>0.1</width_delta><!-- in [fm] -->
       </CausalLiquefier>
   </Liquefier>
-  
+
   <!--Eloss Modules -->
   <Eloss>
       <deltaT>0.1</deltaT>
@@ -92,15 +91,15 @@
       </Lbt>
     <AddLiquefier> true </AddLiquefier>
   </Eloss>
-  
+
   <!-- Hydro  Module 2 -->
   <Hydro>
     <MUSIC>
       <name>MUSIC_2</name>
     </MUSIC>
     <AddLiquefier> true </AddLiquefier>
   </Hydro>
-  
+
   <!-- Jet Hadronization Module -->
   <JetHadronization>
     <name>colorless</name>
@@ -112,5 +111,5 @@
   <SoftParticlization>
     <iSS> </iSS>
   </SoftParticlization>
-  
+
 </jetscape>

external_packages/smash/smash_config.yaml

@@ -2,22 +2,29 @@ Logging:
     default: INFO
 
 General:
-    Modus:  Afterburner   # Dummy, but some value has to be present
+    Modus: List
     Time_Step_Mode: Fixed
-    Delta_Time:    0.1
-    End_Time:      100.0  # Dummy
+    Delta_Time:    0.5
+    End_Time:      100.0  # Reset from the JetScape xml config
     Randomseed:    -1     # Reset from the JetScape xml config
     Nevents:       1      # Dummy
 
-#Collision_Term:
-#    Strings: False
+Modi:
+    List:
+        File_Directory: .
+        File_Prefix: "dummy"
+        Shift_Id: 0
 
-Output:
-    Output_Interval:  5.0
 
-# This section is dummy
-List:
-    File_Directory: .
-    File_Prefix: "dummy"
-    Shift_Id: 0
-    Start_Time: 0.0
+# if some SMASH output in OSCAR,... format is wanted, then uncomment
+# the following lines and add a 'smash_output' folder to your build folder
+Output:
+    Output_Interval:  100.0
+#    Particles:
+#        Format:    ["Oscar2013"]
+#        Extended: True
+#        Only_Final: Yes
+#    Collisions:
+#        Format:    ["Oscar2013"]
+#        Extended: True
+#        Print_Start_End: True

src/initialstate/PythiaGun.cc

@@ -128,6 +128,10 @@ void PythiaGun::InitTask() {
   numbf << eCM;
   readString(numbf.str());
 
+  //Reading vir_factor from xml for MATTER
+  vir_factor = GetXMLElementDouble({"Eloss", "Matter", "vir_factor"});
+  softMomentumCutoff = GetXMLElementDouble({"Hard", "PythiaGun", "softMomentumCutoff"});
+
   std::stringstream lines;
   lines << GetXMLElementText({"Hard", "PythiaGun", "LinesToRead"}, false);
   int i = 0;
@@ -150,8 +154,6 @@ void PythiaGun::InitTask() {
 void PythiaGun::ExecuteTask() {
   VERBOSE(1) << "Run Hard Process : " << GetId() << " ...";
   VERBOSE(8) << "Current Event #" << GetCurrentEvent();
-  //Reading vir_factor from xml for MATTER
-  double vir_factor = GetXMLElementDouble({"Eloss", "Matter", "vir_factor"});
 
   bool flag62 = false;
   vector<Pythia8::Particle> p62;
@@ -197,8 +199,15 @@ void PythiaGun::ExecuteTask() {
 
         // reject rare cases of very soft particles that don't have enough e to get
         // reasonable virtuality
-        if (particle.pT() < 1.0 / sqrt(vir_factor))
+        if (vir_factor > 0. && (particle.pT() < softMomentumCutoff)) {
+          // this cutoff was 1.0/sqrt(vir_factor) in versions < 3.6
+          continue;
+        } else if(vir_factor < 0. && (particle.pAbs() < softMomentumCutoff)) {
           continue;
+        } else if(vir_factor < rounding_error) {
+          JSWARN << "vir_factor should not be zero.";
+          exit(1);
+        }
 
         //if(particle.id()==22) cout<<"########this is a photon!######" <<endl;
         // accept

src/initialstate/PythiaGun.h

@@ -30,6 +30,8 @@ class PythiaGun : public HardProcess, public Pythia8::Pythia {
   double pTHatMin;
   double pTHatMax;
   double eCM;
+  double vir_factor;
+  double softMomentumCutoff;
   bool FSR_on;
   int flag_useHybridHad;
 

src/jet/Matter.cc

@@ -506,7 +506,7 @@ void Matter::DoEnergyLoss(double deltaT, double time, double Q2,
       double max_vir;
       if (vir_factor < 0.0)
         max_vir =
-            pIn[i].e() * pIn[i].e() - pIn[i].restmass() * pIn[i].restmass();
+            std::abs(vir_factor) * (pIn[i].e() * pIn[i].e() - pIn[i].restmass() * pIn[i].restmass());
       else
         max_vir = pT2 * vir_factor;
 
@@ -816,7 +816,7 @@ void Matter::DoEnergyLoss(double deltaT, double time, double Q2,
 		}*/
 
 	    //GeneralQhatFunction(int QhatParametrizationType, double Temperature, double EntropyDensity, double FixAlphas,  double Qhat0, double E, double muSquare)
-	    double muSquare= pIn[i].t(); //Virtuality of the parent; Revisit this when q-hat is virtuality dependent
+	    double muSquare= pIn[i].t(); //Virtuality of the parent; Revist this when q-hat is virtuality dependent
 	    qhatLoc= GeneralQhatFunction(QhatParametrizationType, tempLoc, sdLoc, alphas, qhat0, enerLoc, muSquare);
 
           } else { // outside the QGP medium
@@ -2175,10 +2175,10 @@ double Matter::generate_vac_t_w_M(int p_id, double M, double nu, double t0,
     }
 
     if (std::abs(p_id) == cid || std::abs(p_id) == bid)
-      exit_condition = (std::abs(diff) < s_approx) &&
+      exit_condition = (std::abs(diff) < s_approx) ||
                        (std::abs(t_hi_M0 - t_low_M0) / t_hi_M0 < s_error);
     if (p_id == gid)
-      exit_condition = (std::abs(diff) < s_approx) &&
+      exit_condition = (std::abs(diff) < s_approx) ||
                        (std::abs(t_hi_00 - t_low_00) / t_hi_00 < s_error);
     // need to think about the second statement in the gluon exit condition.
 
@@ -4103,12 +4103,13 @@ double Matter::ModifiedProbability(int QhatParametrization, double tempLoc, doub
       ModifiedAlphas = solve_alphas(qhatLoc, enerLoc, tempLoc);
       break;
 
-      //For HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
-      //Function is 1 / (1+A*pow(log(Q^2),2)+B*pow(log(Q^2),4))
+      //HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
+      //Function is 1/(1+A*pow(log(Q^2),2)+B*pow(log(Q^2),4))
     case 5:
       ModifiedAlphas = RunningAlphaS(ScaleNet)*VirtualityQhatFunction(5,  enerLoc, muSquare) ;
       break;
-      //HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat 
+
+      //HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
       //Function is int^{1}_{xB} e^{-ax} / (1+A*pow(log(Q^2),1)+B*pow(log(Q^2),2))
     case 6:
       ModifiedAlphas = RunningAlphaS(ScaleNet)*VirtualityQhatFunction(6,  enerLoc, muSquare) ;