Main Page   Compound List   File List   Compound Members   File Members  

MyGdDecayModel.cpp

00001 #include <math.h>
00002 #include <iostream.h>
00003 #include "ReactorConstants.hh"
00004 #include "ReactorFortran.hh"
00005 
00006 void MyGdDecayModel(double emax,int& ngamma,double* pgamma){
00007   //
00008   // Ad-hoc model for generating gammas from Gd n-capture.
00009   //
00010   // first pick the number of gammas
00011   //
00012   float xgamma = GammasPerGd-2;
00013   if(xgamma<0)xgamma=0.1;
00014   int error;
00015   rnpssn_(xgamma,ngamma,error);
00016   ngamma+=2;
00017   //
00018   // Assume recoil nucleus can absorn any momentum imbalance and
00019   // just partition the energy among the gammas with directions 
00020   // thrown isotropically.
00021   //  
00022   double* pgptr = pgamma;
00023   float* rv = new float[2*ngamma];
00024   ranlux_(rv,2*ngamma);
00025   for(int i=0;i<ngamma;i++){
00026     double cosg = -1+2*rv[2*i];
00027     double phig = 2*PI*rv[2*i+1];
00028     double t[3];
00029     t[0] = cos(phig)*sqrt(1-cosg*cosg);
00030     t[1] = sin(phig)*sqrt(1-cosg*cosg);
00031     t[2] = cosg;
00032     for(int j=0;j<3;j++)*pgptr++ = t[j];
00033   }
00034   //
00035   // Generate the energies via phase space.
00036   // This involves sampling x(1)*x(2)*..*x(N-1)*(1-x(1)-x(2)-...-x(N-1))
00037   //
00038   double* x = new double[ngamma];
00039   bool OK=0;
00040   while (!OK){
00041     ranlux_(rv,ngamma);
00042     double maxvalue=1;
00043     double logweight=0;
00044     for(int i=0;i<ngamma-1;i++){
00045       x[i]=maxvalue*rv[i];
00046       maxvalue-=x[i];
00047       logweight+=log(x[i]);
00048     }
00049     x[ngamma-1]=maxvalue;
00050     logweight+=log(maxvalue);
00051     double logmaxweight = -(ngamma-1)*log((double) ngamma-1); 
00052     double logtest = log(rv[ngamma-1])+logmaxweight;
00053     OK=(logtest<=logweight);
00054   }
00055   //
00056   // Add energies to gamma 4-vectors
00057   //
00058   pgptr = pgamma;
00059   for(int i=0;i<ngamma;i++){
00060     for(int j=0;j<3;j++){
00061       *pgptr++ *=  emax*x[i];
00062     }
00063   }
00064   delete[] rv;
00065   delete[] x;
00066 }

Generated on Mon Feb 21 16:11:19 2005 for ReactorFsim by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002