In This Package:

Ge68.cc File Reference

#include <stdio.h>
#include <iostream>
#include <CLHEP/Vector/ThreeVector.h>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Units/PhysicalConstants.h>
#include <TRandom.h>
#include <TMath.h>
#include <TF1.h>
#include <TH1.h>

Include dependency graph for Ge68.cc:

Go to the source code of this file.

Namespaces
namespace	std
namespace	CLHEP
Functions
Double_t	dNdE (Double_t *x, Double_t *par)
void	ProcessArgs (int argc, char **argv, long *rseed, char **outfilename, unsigned int *nevents)
void	Usage (char *name)
int	main (int argc, char **argv)
Variables
const double	electronMass = 0.510998910e-3
const double	Q_gs = 1.8991e-3
const double	Q_1st = 0.8218e-3
const double	E_gamma = 1.07735e-3

---

Function Documentation

Double_t dNdE	(	Double_t *	x,
		Double_t *	par
	)

Definition at line 41 of file Ge68.cc.

{
  // par[0] is the Q value in GeV
  // x[0] is electron energy with unit GeV
  Double_t KE = x[0];
  Double_t Q = par[0];

  Double_t Energy = KE + electronMass;
  Double_t p = sqrt(KE*(KE+2*electronMass));

  Double_t spec = (Q-KE)*(Q-KE)*Energy*p;
  return spec;
}

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

Definition at line 157 of file Ge68.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(argv[0]);
      exit(0);
    }
  }
}

void Usage

(

char *

name

)

Definition at line 178 of file Ge68.cc.

                       {
  printf("%s [-seed seed] [-o outputfilename] [-n nevents]\n", name);
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 59 of file Ge68.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);
  
  FILE* stream = stdout;
  if( outFilename!=NULL ) {
    stream = fopen(outFilename, "w");
    if( stream==NULL ) {
      printf("Please enetr a valid filename.\n");
      Usage(argv[0]);
      exit(0);
    }
  }

  gRandom->SetSeed(rseed);

  //create beta spectrum function
  TF1* funcBetaEnergy = new TF1("funcBetaEnergy", dNdE, 0., 2e-3, 1);//GeV


  //HepRandom::setTheSeed(rseed);
  //start by printing some information to comment lines
  fprintf(stream, "# File generated by %s.\n", argv[0]);
  fprintf(stream, "# Random seed for generator = %ld.\n", rseed);
  


  //create a branch weighting table for random sampling
  TH1I *hBranching = new TH1I("hBranching","", 3, 0, 3);
  //See notes at the beginning of the code for the prob values below
  hBranching->SetBinContent(1,87.94);
  hBranching->SetBinContent(2,1.20);
  hBranching->SetBinContent(3,2.02);//lump EC 1st excited state with other cascade
  
  for( size_t i=0 ; i<nEvents ; i++ ) {
    Int_t ibranch = hBranching->FindBin(hBranching->GetRandom());
    double energyInGeV(0.0), momentumInGeV(0.0);
    double px = 0.0, py = 0.0, pz = 0.0;
    double phi = 0.0, theta = 0.0;
    
    if(ibranch==1){//ground state positron emission
      funcBetaEnergy->SetParameter(0, Q_gs);//Q value in GeV for G.S.
      energyInGeV=funcBetaEnergy->GetRandom();
      momentumInGeV = TMath::Sqrt(energyInGeV*energyInGeV
                                  +2.0*electronMass*energyInGeV);      
      phi= gRandom->Uniform(0.0,2.0*TMath::Pi());
      theta= acos(gRandom->Uniform(-1.0,1.0));

      px = cos(phi)*sin(theta)*momentumInGeV;
      py = sin(phi)*sin(theta)*momentumInGeV;
      pz = cos(theta)*momentumInGeV;

      fprintf(stream, "1\n");
      fprintf(stream, "1\t-11 0 0 %e %e %e %e\n", px, py, pz, electronMass); 

    } else if (ibranch==2) {//1st excited state, positron emission
      funcBetaEnergy->SetParameter(0, Q_1st);//Q value in GeV for 1st ES
      energyInGeV=funcBetaEnergy->GetRandom();
      momentumInGeV = TMath::Sqrt(energyInGeV*energyInGeV
                                  +2.0*electronMass*energyInGeV);      
      phi= gRandom->Uniform(0.0,2.0*TMath::Pi());
      theta= acos(gRandom->Uniform(-1.0,1.0));

      px = cos(phi)*sin(theta)*momentumInGeV;
      py = sin(phi)*sin(theta)*momentumInGeV;
      pz = cos(theta)*momentumInGeV;

      fprintf(stream, "2\n");
      fprintf(stream, "1\t-11 0 0 %e %e %e %e\n", px, py, pz, electronMass);      
      //now do the gamma, assuming no angular correlation with e+
      phi= gRandom->Uniform(0.0,2.0*TMath::Pi());
      theta= acos(gRandom->Uniform(-1.0,1.0));
      
      px = cos(phi)*sin(theta)*E_gamma;
      py = sin(phi)*sin(theta)*E_gamma;
      pz = cos(theta)*E_gamma;
      
      fprintf(stream, "1\t22 0 0 %e %e %e 0\n", px, py, pz); 
    } else if (ibranch==3) {//single gamma emission
      phi= gRandom->Uniform(0.0,2.0*TMath::Pi());
      theta= acos(gRandom->Uniform(-1.0,1.0));
      
      px = cos(phi)*sin(theta)*E_gamma;
      py = sin(phi)*sin(theta)*E_gamma;
      pz = cos(theta)*E_gamma;
      fprintf(stream, "1\n");
      fprintf(stream, "1\t22 0 0 %e %e %e 0\n", px, py, pz);
    }
    else {
      cout<<"Wrong branch ID. Ignored ..."<<endl;
    }
  }
  return 0;
}

---

Variable Documentation

const double electronMass = 0.510998910e-3

Definition at line 36 of file Ge68.cc.

const double Q_gs = 1.8991e-3

Definition at line 37 of file Ge68.cc.

const double Q_1st = 0.8218e-3

Definition at line 38 of file Ge68.cc.

const double E_gamma = 1.07735e-3

Definition at line 39 of file Ge68.cc.

---