In This Package:

Ge68.cc

Go to the documentation of this file.

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

using namespace std;
using namespace CLHEP;

//a few hard-coded constants
const double electronMass = 0.510998910e-3; //in GeV 
const double Q_gs = 1.8991e-3; //in GeV
const double Q_1st = 0.8218e-3; //in GeV
const double E_gamma = 1.07735e-3; //in GeV

Double_t dNdE(Double_t *x, Double_t *par)
{
  // 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 );
void Usage(char* name);

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


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(argv[0]);
      exit(0);
    }
  }
}

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

---