Domenico Tiziani / gammalib

Commit 108f91a0697238d1efb7d9bc4eca49df75549b44

Authored by knodlseder
1 parent c266bae6

implement mc() method

Inline Side-by-side
Showing 1 changed file with 30 additions and 8 deletions
src/model/GModelRadialShell.cpp
  View file @108f91a
... ... @@ -40,6 +40,7 @@ const GModelRadialRegistry g_radial_shell_registry(&g_radial_shell_seed);
40 40 /* __ Coding definitions _________________________________________________ */
41 41  
42 42 /* __ Debug definitions __________________________________________________ */
  43 +#define G_DEBUG_MC 0 //!< Debug MC method
43 44  
44 45  
45 46 /*==========================================================================
... ... @@ -316,19 +317,40 @@ double GModelRadialShell::eval_gradients(const double&amp; theta) const
316 317 * @param[in] ran Random number generator.
317 318 *
318 319 * Draws an arbitrary sky position from the 2D shell distribution.
319   - *
320   - * @todo To be implemented.
321 320 ***************************************************************************/
322 321 GSkyDir GModelRadialShell::mc(GRan& ran) const
323 322 {
324   - // Simulate offset from photon arrival direction
325   - // @todo How to draw theta for this model?
326   - double theta = (radius()+width()) * sqrt(ran.uniform());
327   - double phi = 360.0 * ran.uniform();
  323 + // Update precomputation cache
  324 + update();
328 325  
329   - // Rotate sky direction by offset
  326 + // Simulate offset from photon arrival direction
  327 + #if G_DEBUG_MC
  328 + int n_samples = 0;
  329 + #endif
  330 + double theta_max = this->theta_max();
  331 + double sin_theta_max = std::sin(theta_max);
  332 + double u_max = m_norm * sin_theta_max * sin_theta_max;
  333 + double value = 0.0;
  334 + double u = 1.0;
  335 + double theta = 0.0;
  336 + do {
  337 + theta = ran.uniform() * theta_max;
  338 + value = eval(theta) * std::sin(theta);
  339 + u = ran.uniform() * u_max;
  340 + #if G_DEBUG_MC
  341 + n_samples++;
  342 + #endif
  343 + } while (u > value);
  344 + #if G_DEBUG_MC
  345 + std::cout << "#=" << n_samples << " ";
  346 + #endif
  347 +
  348 + // Simulate azimuth angle
  349 + double phi = 360.0 * ran.uniform();
  350 +
  351 + // Rotate pointing direction by offset and azimuth angle
330 352 GSkyDir sky_dir = dir();
331   - sky_dir.rotate(phi, theta);
  353 + sky_dir.rotate(phi, theta*rad2deg);
332 354  
333 355 // Return sky direction
334 356 return sky_dir;
... ...