In This Package:

C11DecayGen.cc

Go to the documentation of this file.

// This is a HepEVT generator of C11 beta+ decay (resulting positrons  ). 
// // Modified from K-40 generator by Raymond Tsang, h0237236@hkusua.hku.hk, 30 Aug, 2006.
// //
// // //\author Fengpeng An
// // //anfp@ihep.ac.cn, 8 Nov., 2009
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <math.h>

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

#include "TFile.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TMath.h"

using namespace std;
using namespace CLHEP;

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


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
  
  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);
   gRandom->SetSeed(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);


// The root file used by this generator should be located at 
// $SITEROOT/dybgaudi/Generators/RadioActivity/C11/
    std::string specfilename = "C11BetaSpectrum.root";
    const char* prefix = getenv("SITEROOT");
    if(prefix) {
        if(strlen(prefix)>0) {
            std::string newname = prefix;
            newname += "/dybgaudi/Generators/RadioActivity/C11/";
            newname += specfilename;
            specfilename = newname;
        }
    }


//cout<<specfilename<<endl;

  TFile  *betaSpec = TFile::Open(specfilename.c_str());

    if(!betaSpec)
        cout << "Can not open C11BetaSpectrum.root. "
               << "Make sure SITEROOT or C11BetaSpectrum.root is set to it's location."
               << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/C11/" << endl;
    if(!(betaSpec->IsOpen()))
        cout << "Can not open C11BetaSpectrum.root "
               << "Make sure SITEROOT is set to it's location."
               << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/C11/" << endl;



  double cosTheta, theta, azimuth; // angles used
  unsigned int i;

  TH1F *h1 = (TH1F*)betaSpec->Get("C11");
  
  for( i=0 ; i<nEvents ; i++ ) {

      cosTheta = RandFlat::shoot(-1, 1);
      azimuth = RandFlat::shoot( 2.0*M_PI );
      theta = acos(cosTheta);
      p1.setRThetaPhi(1, acos(cosTheta), azimuth);

      double ke = h1->GetRandom();
      double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
      p1=p1*pe*keV;
      fprintf(stream, "1\n");
     fprintf(stream, "1\t-11 0 0 %e %e %e %e\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV, electron_mass_c2/GeV); 

  
  }

  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("C11.exe [-seed seed] [-o outputfilename] [-n nevents]\n");
}

---