---

MyGdDecayModel.cpp

#include <math.h>
#include <iostream.h>
#include "ReactorConstants.hh"
#include "ReactorFortran.hh"

void MyGdDecayModel(double emax,int& ngamma,double* pgamma){
  //
  // Ad-hoc model for generating gammas from Gd n-capture.
  //
  static ReactorConstants k;
  double pi = k.pi;
  //
  // first pick the number of gammas
  //
  float xgamma = k.GammasPerGd-2;
  if(xgamma<0)xgamma=0.1;
  int error;
  rnpssn_(xgamma,ngamma,error);
  ngamma+=2;
  //
  // Assume recoil nucleus can absorn any momentum imbalance and
  // just partition the energy among the gammas with directions 
  // thrown isotropically.
  //  
  double* pgptr = pgamma;
  float rv[2*ngamma];
  ranlux_(rv,2*ngamma);
  for(int i=0;i<ngamma;i++){
    double cosg = -1+2*rv[2*i];
    double phig = 2*pi*rv[2*i+1];
    double t[3];
    t[0] = cos(phig)*sqrt(1-cosg*cosg);
    t[1] = sin(phig)*sqrt(1-cosg*cosg);
    t[2] = cosg;
    for(int j=0;j<3;j++)*pgptr++ = t[j];
  }
  //
  // Generate the energies via phase space.
  // This involves sampling x(1)*x(2)*..*x(N-1)*(1-x(1)-x(2)-...-x(N-1))
  //
  double x[ngamma];
  bool OK=0;
  while (!OK){
    ranlux_(rv,ngamma);
    double maxvalue=1;
    double logweight=0;
    for(int i=0;i<ngamma-1;i++){
      x[i]=maxvalue*rv[i];
      maxvalue-=x[i];
      logweight+=log(x[i]);
    }
    x[ngamma-1]=maxvalue;
    logweight+=log(maxvalue);
    double logmaxweight = -(ngamma-1)*log(ngamma-1); 
    double logtest = log(rv[ngamma-1])+logmaxweight;
    OK=(logtest<=logweight);
  }
  //
  // Add energies to gamma 4-vectors
  //
  pgptr = pgamma;
  for(int i=0;i<ngamma;i++){
    for(int j=0;j<3;j++){
      *pgptr++ *=  emax*x[i];
    }
  }
}

---