/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Potassium_40_gammas.cc
Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 #include <cstdio>
00010 #include <cstring>
00011 #include <math.h>
00012 #include <iostream>
00013 #include "TMath.h"
00014 
00015 #include <CLHEP/Vector/ThreeVector.h>
00016 #include <CLHEP/Random/Randomize.h>
00017 #include <CLHEP/Units/PhysicalConstants.h>
00018 
00019 
00020 using namespace std;
00021 using namespace CLHEP;
00022 
00023 void ProcessArgs(int argc, char** argv, long* rseed, char** outfilename, 
00024                  unsigned int* nevents );
00025 void Usage();
00026 
00027 //------------------------------------------------------------
00028 //beta spectrum is from GEANT4
00029 class BetaDecay {
00030     int A, Z;
00031     double Q, T_endpoint, PI;
00032 
00033     double norm;
00034 
00035     bool doff, onlyff;
00036 
00037 public:
00038     // Q is unitless decay energy, t_endpoint is measured unitless
00039     // kinetic energy endpoint of beta.
00040     BetaDecay(int a, int z, double q, double t_endpoint) 
00041         : A(a), Z(z), Q(q), T_endpoint(t_endpoint), PI(TMath::Pi()), norm(0)
00042     {
00043         doff = true;            // for testing w/out FF
00044         onlyff = false;         // for testing only FF
00045         norm = this->dnde_norm(t_endpoint);
00046     }
00047 
00048     // Normalize dnde_noff (find peak)
00049     double dnde_norm(double T_end) {
00050         const int nbins = 100;
00051         const double dT = T_end/nbins;
00052         double sum = 0.0;
00053         for (int ind = 1; ind <= nbins; ++ind) {
00054             double T = double(ind)*dT;
00055             double dnde = 0.0;
00056             if (onlyff) dnde = this->GetFF(T);
00057             else{
00058                 dnde = this->dnde_noff(T);
00059                 if (doff) dnde *= this->GetFF(T);
00060             }
00061             sum += dnde;
00062         }
00063         return sum/nbins;
00064     }
00065 
00066     double dnde(double T) {
00067         double dnde = 0.0;
00068         if (onlyff) dnde = this->GetFF(T);
00069         else {
00070             dnde = this->dnde_noff(T);
00071             if (doff) dnde *= this->GetFF(T);
00072         }
00073         return  dnde / norm;
00074     }
00075 
00076     // Beta spectrum ignoring nuclear coulomb field. T = unitless
00077     // kinetic energy.
00078     double dnde_noff(double T) {
00079         double ee = Q+1;
00080         //double ee = T_endpoint+1;
00081         double g = T+1.;
00082         return std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
00083     }
00084 
00085     // Fermi function to correct for nuclear coulomb field. T =
00086     // unitless kinetic energy of beta.  Unchanged except to replace
00087     // "G4" types with std ones and rename input var.
00088     double GetFF(double T) {
00089         double A1, A2;
00090         double P, U, S, Y;
00091         double F2;
00092         double E = T+1.;  
00093         P=std::sqrt(E*E-1.0) ;
00094         U=Z/137.0;
00095         S=std::sqrt(1.0-U*U) - 1.;
00096         Y = 2*PI*U*E/P;
00097         A1 = U*U*E*E + P*P/4.;
00098         A2 = std::fabs(Y/(1-std::exp(-Y)));
00099         F2 = std::pow(A1,S) * A2; 
00100         return F2;
00101     }
00102 
00103     // Normalize spectrum from GetFF (find peak).  T_end = unitless
00104     // end point kinetic energy of beta.  unchanged except to replace
00105     // "G4" types with std ones and rename input var.
00106     double GetFFN(double T_end) {
00107         double A1, A2;
00108         double P, U, S, Y;
00109         double F2,E;
00110         double EE = T_end/100.;
00111         U=Z/137.0;
00112         S=std::sqrt(1.0-U*U) - 1.;
00113         double F1 = 1E-10;
00114         for (int i = 1; i<=100 ; i++) {
00115             E = double(i)*EE + 1.;
00116             P=std::sqrt(E*E-1.0) ;
00117             Y = 2*PI*U*E/P;
00118             A1 = U*U*E*E + P*P/4.;
00119             A2 = std::fabs(Y/(1-std::exp(-Y)));
00120             F2 = std::pow(A1,S) * A2; 
00121             if (F2 > F1) F1 = F2;
00122         }                  
00123         return F1;
00124     }
00125 
00126 };
00127 //-------------------------------------------------------------------
00128 
00129 int main(int argc, char** argv) {
00130   long rseed = 0;
00131   char* outFilename = NULL;
00132   unsigned int nEvents = 1000000000; // a billion. Default to something big for piping 
00133   ProcessArgs(argc, argv, &rseed, &outFilename, &nEvents);
00134   Hep3Vector p1; // the gamma momentum vector. Unit GeV/c
00135   Hep3Vector p2; // the electron momentum vector. Unit GeV/c
00136   
00137   FILE* stream = stdout;
00138   if( outFilename!=NULL ) {
00139     stream = fopen(outFilename, "w");
00140     if( stream==NULL ) {
00141       printf("Please enter a valid filename.\n");
00142       Usage();
00143       exit(0);
00144     }
00145   }
00146   HepRandom::setTheSeed(rseed);
00147   //start by printing some information to comment lines
00148   fprintf(stream, "# File generated by %s.\n", argv[0]);
00149   fprintf(stream, "# Ransom seed for generator = %ld.\n", rseed);
00150   
00151   double cosTheta, theta, azimuth; // angles used
00152   unsigned int i;
00153   double bf; //K40 decay branching fraction
00154   for( i=0 ; i<nEvents ; i++ ) {
00155     cosTheta = RandFlat::shoot(-1., 1.);
00156     azimuth = RandFlat::shoot( 2.0*M_PI );
00157     theta = acos(cosTheta);
00158 
00159     bf = RandFlat::shoot(0.,1.);
00160     if(bf<0.107){  //10.7% gamma radiation
00161       p1.setRThetaPhi(1, theta, azimuth);
00162       p1*=1.461*MeV;
00163       fprintf(stream, "1\n");
00164       fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV );
00165     }
00166     else{ //89.3% beta decay, maximum energy is 1.31 MeV.
00167       bool accept = 0;
00168       double emass = 0.51099892;
00169       while(!accept){
00170         double rand1 = RandFlat::shoot(0.,1.31);//random energy
00171         double rand2 = RandFlat::shoot(0.,1.7);//1.7> spectrum max
00172         BetaDecay bd = BetaDecay(0,+20,1.31/emass,1.31/emass);
00173         if(bd.dnde(rand1/emass)>rand2) {//accept
00174           //rand1 is kinetic energy, now translate into momentum
00175           double momentum = sqrt((rand1+emass)*(rand1+emass)-emass*emass);
00176           p2.setRThetaPhi(1, theta, azimuth);
00177           p2*=momentum*MeV;
00178 
00179           //electron PDG id is 11
00180           fprintf(stream, "1\n"); 
00181           fprintf(stream, "1\t11 0 0 %e %e %e 0.51099892e-03\n", p2.x()/GeV, p2.y()/GeV, p2.z()/GeV );
00182 
00183           accept=1;
00184         }
00185       }   
00186     }
00187   }
00188   return 0;
00189 }
00190 
00191 
00192 void ProcessArgs(int argc, char** argv, long *rseed, 
00193                  char** outfilename, unsigned int* nevents) {
00194   int i;
00195   for( i=1 ; i<argc ; i++ ) {
00196     if( strcmp(argv[i], "-seed")==0 ) {
00197       i++;
00198       sscanf(argv[i], "%ld", rseed);
00199     } else if( strcmp(argv[i], "-o")==0 ) {
00200       i++;
00201       *outfilename = new char[strlen(argv[i]) +1];
00202       strcpy(*outfilename, argv[i]);
00203     } else if( strcmp(argv[i], "-n")==0 ) { 
00204       i++;
00205       sscanf(argv[i], "%ud", nevents);
00206     } else {
00207       Usage();
00208       exit(0);
00209     }
00210   }
00211 }
00212 
00213 void Usage() {
00214   printf("Potassium-40 [-seed seed] [-o outputfilename] [-n nevents]\n");
00215 }
00216 
00217 
00218 
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 09:52:24 for K40 by doxygen 1.7.4