In This Package:

Potassium_40_gammas.cc File Reference

#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>

Include dependency graph for Potassium_40_gammas.cc:

Go to the source code of this file.

Namespaces
namespace	std
namespace	CLHEP
Classes
class	BetaDecay
Functions
void	ProcessArgs (int argc, char **argv, long *rseed, char **outfilename, unsigned int *nevents)
void	Usage ()
int	main (int argc, char **argv)

---

Function Documentation

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

Definition at line 192 of file Potassium_40_gammas.cc.

                                                            {
  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

(

)

Definition at line 213 of file Potassium_40_gammas.cc.

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

int main	(	int	argc,
		char **	argv
	)

Definition at line 129 of file Potassium_40_gammas.cc.

                                {
  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;
}

---