trying to get pdf from inelastic disrtibutions

openmc-dev · itay-space · Oct 25, 2022 · Oct 25, 2022 · Oct 25, 2022 · Oct 25, 2022
commit 60df64378fe8999991e1f5c2984e0fd7b1fae807
diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h
@@ -43,7 +43,7 @@ class ReactionProduct {
   //! \param[out] mu Outgoing cosine with respect to current direction
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
-
+  double get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const;
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
   double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]

diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h
@@ -40,6 +40,8 @@ class CorrelatedAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+ void get_pdf(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
   // energy property
   vector<double>& energy() { return energy_; }

diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h
@@ -31,7 +31,7 @@ class KalbachMann : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
-   double get_pdf(double E_in,double mymu, uint64_t* seed) const;
+   double get_pdf(double E_in,double& E_out,double mymu, uint64_t* seed) const;
 
 private:
   //! Outgoing energy/angle at a single incoming energy

diff --git a/src/physics.cpp b/src/physics.cpp
@@ -1110,6 +1110,7 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
   // sample outgoing energy and scattering cosine
   double E;
   double mu;
+  //std::cout << "E_in inelastic_scatter" << E_in << std::endl;
   rx.products_[0].sample(E_in, E, mu, p.current_seed());
 
   // if scattering system is in center-of-mass, transfer cosine of scattering

diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
@@ -108,4 +108,77 @@ void ReactionProduct::sample(
   }
 }
 
+double ReactionProduct::get_pdf(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  double mypdf = 0;
+
+  int distribution_index;
+  auto n = applicability_.size();
+  if (n > 1) {
+    double prob = 0.0;
+    double c = prn(seed);
+    for (int i = 0; i < n; ++i) {
+      // Determine probability that i-th energy distribution is sampled
+      prob += applicability_[i](E_in);
+
+      // If i-th distribution is sampled, sample energy from the distribution
+      if (c <= prob) {
+        //distribution_[i]->sample(E_in, E_out, mu, seed);
+        distribution_index = i;
+        break;
+      }
+    }
+  } else {
+    // If only one distribution is present, go ahead and sample it
+    //distribution_[0]->sample(E_in, E_out, mu, seed);
+    distribution_index = 0;
+
+  }
+ // now extract pdf 
+
+AngleEnergy* angleEnergyPtr = distribution_[distribution_index].get();
+
+if (CorrelatedAngleEnergy* correlatedAE = dynamic_cast<CorrelatedAngleEnergy*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*correlatedAE).name() << " implementation." << std::endl;
+    // Handle CorrelatedAngleEnergy
+} else if (KalbachMann* kalbachMann = dynamic_cast<KalbachMann*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*kalbachMann).name() << " implementation." << std::endl;
+    mypdf = (*kalbachMann).get_pdf(E_in,E_out, 0.3333 , seed);
+   // std::cout << "mypdf " << (*kalbachMann).get_pdf(E_in,E_out, 0.3333 , seed) << std::endl;
+   // std::cout << " my E_in " << E_in <<std::endl;
+   // std::cout << " my E out " << E_out <<std::endl;
+    // Handle KalbachMann
+} else if (NBodyPhaseSpace* nBodyPS = dynamic_cast<NBodyPhaseSpace*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*nBodyPS).name() << " implementation." << std::endl;
+    // Handle NBodyPhaseSpace
+} else if (UncorrelatedAngleEnergy* uncorrelatedAE = dynamic_cast<UncorrelatedAngleEnergy*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*uncorrelatedAE).name() << " implementation." << std::endl;
+    // Handle UncorrelatedAngleEnergy
+} else if (CoherentElasticAE* coherentElasticAE = dynamic_cast<CoherentElasticAE*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*coherentElasticAE).name() << " implementation." << std::endl;
+    // Handle CoherentElasticAE
+} else if (IncoherentElasticAE* incoherentElasticAE = dynamic_cast<IncoherentElasticAE*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*incoherentElasticAE).name() << " implementation." << std::endl;
+    // Handle IncoherentElasticAE
+} else if (IncoherentElasticAEDiscrete* incoherentElasticAEDiscrete = dynamic_cast<IncoherentElasticAEDiscrete*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*incoherentElasticAEDiscrete).name() << " implementation." << std::endl;
+    // Handle IncoherentElasticAEDiscrete
+} else if (IncoherentInelasticAEDiscrete* incoherentInelasticAEDiscrete = dynamic_cast<IncoherentInelasticAEDiscrete*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*incoherentInelasticAEDiscrete).name() << " implementation." << std::endl;
+    // Handle IncoherentInelasticAEDiscrete
+} else if (IncoherentInelasticAE* incoherentInelasticAE = dynamic_cast<IncoherentInelasticAE*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*incoherentInelasticAE).name() << " implementation." << std::endl;
+    // Handle IncoherentInelasticAE
+} else if (MixedElasticAE* mixedElasticAE = dynamic_cast<MixedElasticAE*>(angleEnergyPtr)) {
+    std::cout << "Used " << typeid(*mixedElasticAE).name() << " implementation." << std::endl;
+    // Handle MixedElasticAE
+} else {
+    std::cout << "Unknown derived type." << std::endl;
+}
+
+
+ return mypdf;
+}
+
 } // namespace openmc
diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp
@@ -263,4 +263,114 @@ void CorrelatedAngleEnergy::sample(
   }
 }
 
+void CorrelatedAngleEnergy::get_pdf(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  // Find energy bin and calculate interpolation factor -- if the energy is
+  // outside the range of the tabulated energies, choose the first or last bins
+  auto n_energy_in = energy_.size();
+  int i;
+  double r;
+  if (E_in < energy_[0]) {
+    i = 0;
+    r = 0.0;
+  } else if (E_in > energy_[n_energy_in - 1]) {
+    i = n_energy_in - 2;
+    r = 1.0;
+  } else {
+    i = lower_bound_index(energy_.begin(), energy_.end(), E_in);
+    r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);
+  }
+
+  // Sample between the ith and [i+1]th bin
+  int l = r > prn(seed) ? i + 1 : i;
+
+  // Interpolation for energy E1 and EK
+  int n_energy_out = distribution_[i].e_out.size();
+  int n_discrete = distribution_[i].n_discrete;
+  double E_i_1 = distribution_[i].e_out[n_discrete];
+  double E_i_K = distribution_[i].e_out[n_energy_out - 1];
+
+  n_energy_out = distribution_[i + 1].e_out.size();
+  n_discrete = distribution_[i + 1].n_discrete;
+  double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
+  double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
+
+  double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
+  double E_K = E_i_K + r * (E_i1_K - E_i_K);
+
+  // Determine outgoing energy bin
+  n_energy_out = distribution_[l].e_out.size();
+  n_discrete = distribution_[l].n_discrete;
+  double r1 = prn(seed);
+  double c_k = distribution_[l].c[0];
+  int k = 0;
+  int end = n_energy_out - 2;
+
+  // Discrete portion
+  for (int j = 0; j < n_discrete; ++j) {
+    k = j;
+    c_k = distribution_[l].c[k];
+    if (r1 < c_k) {
+      end = j;
+      break;
+    }
+  }
+
+  // Continuous portion
+  double c_k1;
+  for (int j = n_discrete; j < end; ++j) {
+    k = j;
+    c_k1 = distribution_[l].c[k + 1];
+    if (r1 < c_k1)
+      break;
+    k = j + 1;
+    c_k = c_k1;
+  }
+
+  double E_l_k = distribution_[l].e_out[k];
+  double p_l_k = distribution_[l].p[k];
+  if (distribution_[l].interpolation == Interpolation::histogram) {
+    // Histogram interpolation
+    if (p_l_k > 0.0 && k >= n_discrete) {
+      E_out = E_l_k + (r1 - c_k) / p_l_k;
+    } else {
+      E_out = E_l_k;
+    }
+
+  } else if (distribution_[l].interpolation == Interpolation::lin_lin) {
+    // Linear-linear interpolation
+    double E_l_k1 = distribution_[l].e_out[k + 1];
+    double p_l_k1 = distribution_[l].p[k + 1];
+
+    double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
+    if (frac == 0.0) {
+      E_out = E_l_k + (r1 - c_k) / p_l_k;
+    } else {
+      E_out =
+        E_l_k +
+        (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
+          p_l_k) /
+          frac;
+    }
+  }
+
+  // Now interpolate between incident energy bins i and i + 1
+  if (k >= n_discrete) {
+    if (l == i) {
+      E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
+    } else {
+      E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
+    }
+  }
+
+  // Find correlated angular distribution for closest outgoing energy bin
+  if (r1 - c_k < c_k1 - r1 ||
+      distribution_[l].interpolation == Interpolation::histogram) {
+    mu = distribution_[l].angle[k]->sample(seed);
+  } else {
+    mu = distribution_[l].angle[k + 1]->sample(seed);
+  }
+}
+
 } // namespace openmc
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
@@ -237,12 +237,12 @@ void KalbachMann::sample(
     mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
   }
 }
-double KalbachMann::get_pdf(double E_in,double mymu, uint64_t* seed) const
+double KalbachMann::get_pdf(double E_in,double& E_out,double mymu, uint64_t* seed) const
 {
   // Find energy bin and calculate interpolation factor -- if the energy is
   // outside the range of the tabulated energies, choose the first or last bins
   double myE_in = E_in;
-  double myE_out;
+  //double E_out;
   auto n_energy_in = energy_.size();
   int i;
   double r;
@@ -309,14 +309,15 @@ double KalbachMann::get_pdf(double E_in,double mymu, uint64_t* seed) const
   if (distribution_[l].interpolation == Interpolation::histogram) {
     // Histogram interpolation
     if (p_l_k > 0.0 && k >= n_discrete) {
-      myE_out = E_l_k + (r1 - c_k) / p_l_k;
+      E_out = E_l_k + (r1 - c_k) / p_l_k;
     } else {
-      myE_out = E_l_k;
+      E_out = E_l_k;
     }
-
+     std::cout << " Histogram interpolation" <<std::endl;
     // Determine Kalbach-Mann parameters
     km_r = distribution_[l].r[k];
     km_a = distribution_[l].a[k];
+    //std::cout <<"Histo " << " km_r " <<  km_r << " km_a " <<  km_a << std::endl;
 
   } else {
     // Linear-linear interpolation
@@ -325,24 +326,35 @@ double KalbachMann::get_pdf(double E_in,double mymu, uint64_t* seed) const
 
     double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
     if (frac == 0.0) {
-      myE_out = E_l_k + (r1 - c_k) / p_l_k;
+      E_out = E_l_k + (r1 - c_k) / p_l_k;
     } else {
-      myE_out =
+      E_out =
         E_l_k +
         (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
           p_l_k) /
           frac;
     }
-
+   std::cout << " Linear-linear interpolation " <<std::endl;
     // Determine Kalbach-Mann parameters
     km_r = distribution_[l].r[k] +
-           (myE_out - E_l_k) / (E_l_k1 - E_l_k) *
+           (E_out - E_l_k) / (E_l_k1 - E_l_k) *
              (distribution_[l].r[k + 1] - distribution_[l].r[k]);
     km_a = distribution_[l].a[k] +
-           (myE_out - E_l_k) / (E_l_k1 - E_l_k) *
+           (E_out - E_l_k) / (E_l_k1 - E_l_k) *
              (distribution_[l].a[k + 1] - distribution_[l].a[k]);
   }
 
+
+// Now interpolate between incident energy bins i and i + 1
+  if (k >= n_discrete) {
+    if (l == i) {
+      E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
+    } else {
+      E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
+    }
+  }
+  //std::cout << " my E out in KM  " << E_out <<std::endl;
+
   // https://docs.openmc.org/en/v0.8.0/methods/physics.html#equation-KM-pdf-angle
    double pdf_mu = km_a / (2 * std::sinh(km_a)) * (std::cosh(km_a * mymu) + km_r * std::sinh(km_a * mymu));
    return pdf_mu;

diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
@@ -76,7 +76,7 @@ double UncorrelatedAngleEnergy::get_pdf(
     // no angle distribution given => assume isotropic for all energies
     pdf = 0.5;
   }
-
+//E_out = energy_->sample(E_in, seed);
 return pdf;
 }
 

diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
@@ -2717,32 +2717,33 @@ void score_point_tally(Particle& p)
 
 
  if (p.event_mt() != 2){
-   //std::cout << "my mt  "<< p.event_mt() <<std::endl;
-   /*
-   const auto& rx {nuc->reactions_[ p.event_index_mt()]};
-   auto n = applicability_.size();
-   int distribution_index;
-    if (n > 1) {
-    double prob = 0.0;
-    double c = prn(seed);
-    for (int i = 0; i < n; ++i) {
-      // Determine probability that i-th energy distribution is sampled
-      prob += applicability_[i](E_in);
-
-      // If i-th distribution is sampled, sample energy from the distribution
-      if (c <= prob) {
-        distribution_[i]->sample(E_in, E_out, mu, seed)
-        distribution_index = i;
-        break;
-      }
-    }
-  } else {
-    // If only one distribution is present, go ahead and sample it
-    distribution_[0]->sample(E_in, E_out, mu, seed);
-    distribution_index = 0;
+   std::cout << "my mt  "<< p.event_mt() <<std::endl;
+   // kalbach
+   // get E_out from get_pdf
+  //mu_cm -> from E_out and  Jacobian
+  // make sure v_t is 0
+
+//
+
+    // copy energy of neutron
+   double E_in = p.E_last();
+   std::cout << "E_in collision " << E_in<<std::endl;
+   // sample outgoing energy and scattering cosine
+  double E_out;
+  double mu;
+
+
+ const auto& rx {nuc->reactions_[ p.event_index_mt()]};
+ double d = rx->products_[0].get_pdf(E_in, E_out, mu, p.current_seed());
+   if (rx->scatter_in_cm_) 
+  {
+    std::cout << "pdf com " << d <<std::endl;
+    std::cout << "E_out com " << E_out <<std::endl;
   }
+
+// Now check which distribution is used 
+
 
-  */
  }
 
  // starting scroing loop on ghost particles