added v_t = velcity of target particle as a property of particle data

openmc-dev · Oct 25, 2022 · 65f3033 · Itay-max · Oct 26, 2022 · 65f3033
1 parent 913273d
commit 65f3033
Show file tree

Hide file tree

Showing 3 changed files with 9 additions and 2 deletions.
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
@@ -251,6 +251,7 @@ class ParticleData {
   int event_nuclide_;     //!< index in nuclides array
   int event_mt_;          //!< reaction MT
   int delayed_group_ {0}; //!< delayed group
+  Direction v_t_;
 
   // Post-collision physical data
   int n_bank_ {0};        //!< number of fission sites banked
@@ -382,6 +383,8 @@ class ParticleData {
   const int& event_nuclide() const { return event_nuclide_; }
   int& event_mt() { return event_mt_; }
   int& delayed_group() { return delayed_group_; }
+  Position& v_t() { return v_t_; }
+  const Position& v_t() const { return v_t_; }
 
   int& n_bank() { return n_bank_; }
   int& n_bank_second() { return n_bank_second_; }

diff --git a/src/physics.cpp b/src/physics.cpp
@@ -744,7 +744,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p)
     v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n,
       p.neutron_xs(i_nuclide).elastic, kT, p.current_seed());
   }
-
+  p.v_t() = C_LIGHT*std::sqrt(2/(awr*p.getMass()))*v_t;
   // Velocity of center-of-mass
   Direction v_cm = (v_n + awr * v_t) / (awr + 1.0);
 

diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
@@ -2436,7 +2436,7 @@ void score_point_tally(Particle& p)
   col_counter ++;
   fmt::print("------------------------collison happened------------------------\n");
   fmt::print("col counter = {}\n",col_counter);
-  //std::cout << "mass in ev  " << p.getMass() << std::endl ;
+  std::cout << "mass in ev  " << p.getMass() << std::endl ;
   // Determine the collision estimate of the flux
   bool verbose=false;//true;
   double ReturnArray[4]= {std::nan(""),std::nan(""),std::nan(""),std::nan("")};
@@ -3084,6 +3084,10 @@ void getMu_COM(double x_det , double y_det , double z_det ,Particle p_col , doub
   double m3= m1; // mass of outgoing particle to detector  (rest mass?)
   double m4= m2; // mass of recoil target  system
   double p1[3]={p1_tot*p_col.u_last().x,p1_tot* p_col.u_last().y,p1_tot* p_col.u_last().z}; // 3 momentum of incoming particle
+  //p.v_t() * m2 / C_LIGHT;
+  std::cout << "p1:  "<< p1[0] <<" "<< p1[1]<<" " <<p1[2] <<std::endl;
+  std::cout << "p2:  "<< (p_col.v_t() * m2 / C_LIGHT).x <<" " <<(p_col.v_t() * m2 / C_LIGHT).y << " "<<(p_col.v_t() * m2 / C_LIGHT).z <<std::endl;
+  //std::cout << (p.v_t() * m2 / C_LIGHT).x << (p.v_t() * m2 / C_LIGHT).y << (p.v_t() * m2 / C_LIGHT).z <<std::endl;
   double p2[3]={0, 0, 0}; //3 momentum of target in lab  
 
   // calculate