Potassium_40_gammas.cc

Go to the documentation of this file.

#include <cstdio>
#include <cstring>
#include <math.h>
#include <iostream>
#include "TMath.h"

#include <CLHEP/Vector/ThreeVector.h>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Units/PhysicalConstants.h>


using namespace std;
using namespace CLHEP;

void ProcessArgs(int argc, char** argv, long* rseed, char** outfilename, 
                 unsigned int* nevents );
void Usage();

//------------------------------------------------------------
//beta spectrum is from GEANT4
class BetaDecay {
    int A, Z;
    double Q, T_endpoint, PI;

    double norm;

    bool doff, onlyff;

public:
    // Q is unitless decay energy, t_endpoint is measured unitless
    // kinetic energy endpoint of beta.
    BetaDecay(int a, int z, double q, double t_endpoint) 
        : A(a), Z(z), Q(q), T_endpoint(t_endpoint), PI(TMath::Pi()), norm(0)
    {
        doff = true;            // for testing w/out FF
        onlyff = false;         // for testing only FF
        norm = this->dnde_norm(t_endpoint);
    }

    // Normalize dnde_noff (find peak)
    double dnde_norm(double T_end) {
        const int nbins = 100;
        const double dT = T_end/nbins;
        double sum = 0.0;
        for (int ind = 1; ind <= nbins; ++ind) {
            double T = double(ind)*dT;
            double dnde = 0.0;
            if (onlyff) dnde = this->GetFF(T);
            else{
                dnde = this->dnde_noff(T);
                if (doff) dnde *= this->GetFF(T);
            }
            sum += dnde;
        }
        return sum/nbins;
    }

    double dnde(double T) {
        double dnde = 0.0;
        if (onlyff) dnde = this->GetFF(T);
        else {
            dnde = this->dnde_noff(T);
            if (doff) dnde *= this->GetFF(T);
        }
        return  dnde / norm;
    }

    // Beta spectrum ignoring nuclear coulomb field. T = unitless
    // kinetic energy.
    double dnde_noff(double T) {
        double ee = Q+1;
        //double ee = T_endpoint+1;
        double g = T+1.;
        return std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
    }

    // Fermi function to correct for nuclear coulomb field. T =
    // unitless kinetic energy of beta.  Unchanged except to replace
    // "G4" types with std ones and rename input var.
    double GetFF(double T) {
        double A1, A2;
        double P, U, S, Y;
        double F2;
        double E = T+1.;  
        P=std::sqrt(E*E-1.0) ;
        U=Z/137.0;
        S=std::sqrt(1.0-U*U) - 1.;
        Y = 2*PI*U*E/P;
        A1 = U*U*E*E + P*P/4.;
        A2 = std::fabs(Y/(1-std::exp(-Y)));
        F2 = std::pow(A1,S) * A2; 
        return F2;
    }

    // Normalize spectrum from GetFF (find peak).  T_end = unitless
    // end point kinetic energy of beta.  unchanged except to replace
    // "G4" types with std ones and rename input var.
    double GetFFN(double T_end) {
        double A1, A2;
        double P, U, S, Y;
        double F2,E;
        double EE = T_end/100.;
        U=Z/137.0;
        S=std::sqrt(1.0-U*U) - 1.;
        double F1 = 1E-10;
        for (int i = 1; i<=100 ; i++) {
            E = double(i)*EE + 1.;
            P=std::sqrt(E*E-1.0) ;
            Y = 2*PI*U*E/P;
            A1 = U*U*E*E + P*P/4.;
            A2 = std::fabs(Y/(1-std::exp(-Y)));
            F2 = std::pow(A1,S) * A2; 
            if (F2 > F1) F1 = F2;
        }                  
        return F1;
    }

};
//-------------------------------------------------------------------

int main(int argc, char** argv) {
  long rseed = 0;
  char* outFilename = NULL;
  unsigned int nEvents = 1000000000; // a billion. Default to something big for piping 
  ProcessArgs(argc, argv, &rseed, &outFilename, &nEvents);
  Hep3Vector p1; // the gamma momentum vector. Unit GeV/c
  Hep3Vector p2; // the electron momentum vector. Unit GeV/c
  
  FILE* stream = stdout;
  if( outFilename!=NULL ) {
    stream = fopen(outFilename, "w");
    if( stream==NULL ) {
      printf("Please enter a valid filename.\n");
      Usage();
      exit(0);
    }
  }
  HepRandom::setTheSeed(rseed);
  //start by printing some information to comment lines
  fprintf(stream, "# File generated by %s.\n", argv[0]);
  fprintf(stream, "# Ransom seed for generator = %ld.\n", rseed);
  
  double cosTheta, theta, azimuth; // angles used
  unsigned int i;
  double bf; //K40 decay branching fraction
  for( i=0 ; i<nEvents ; i++ ) {
    cosTheta = RandFlat::shoot(-1., 1.);
    azimuth = RandFlat::shoot( 2.0*M_PI );
    theta = acos(cosTheta);

    bf = RandFlat::shoot(0.,1.);
    if(bf<0.107){  //10.7% gamma radiation
      p1.setRThetaPhi(1, theta, azimuth);
      p1*=1.461*MeV;
      fprintf(stream, "1\n");
      fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV );
    }
    else{ //89.3% beta decay, maximum energy is 1.31 MeV.
      bool accept = 0;
      double emass = 0.51099892;
      while(!accept){
        double rand1 = RandFlat::shoot(0.,1.31);//random energy
        double rand2 = RandFlat::shoot(0.,1.7);//1.7> spectrum max
        BetaDecay bd = BetaDecay(0,+20,1.31/emass,1.31/emass);
        if(bd.dnde(rand1/emass)>rand2) {//accept
          //rand1 is kinetic energy, now translate into momentum
          double momentum = sqrt((rand1+emass)*(rand1+emass)-emass*emass);
          p2.setRThetaPhi(1, theta, azimuth);
          p2*=momentum*MeV;

          //electron PDG id is 11
          fprintf(stream, "1\n"); 
          fprintf(stream, "1\t11 0 0 %e %e %e 0.51099892e-03\n", p2.x()/GeV, p2.y()/GeV, p2.z()/GeV );

          accept=1;
        }
      }   
    }
  }
  return 0;
}


void ProcessArgs(int argc, char** argv, long *rseed, 
                 char** outfilename, unsigned int* nevents) {
  int i;
  for( i=1 ; i<argc ; i++ ) {
    if( strcmp(argv[i], "-seed")==0 ) {
      i++;
      sscanf(argv[i], "%ld", rseed);
    } else if( strcmp(argv[i], "-o")==0 ) {
      i++;
      *outfilename = new char[strlen(argv[i]) +1];
      strcpy(*outfilename, argv[i]);
    } else if( strcmp(argv[i], "-n")==0 ) { 
      i++;
      sscanf(argv[i], "%ud", nevents);
    } else {
      Usage();
      exit(0);
    }
  }
}

void Usage() {
  printf("Potassium-40 [-seed seed] [-o outputfilename] [-n nevents]\n");
}