In This Package:

C9DecayGen.cc

Go to the documentation of this file.

// This is a HepEVT generator of C9 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

// data: http://www.nndc.bnl.gov/nudat2/getdataset.jsp?nucleus=9B&unc=nds
#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 p0; // the beta  momentum vector. Unit GeV/c

// some more detailed processes will be added later when necessary. 
  
  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, "# Ransom seed for generator = %ld.\n", rseed);
  

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


//cout<<specfilename<<endl;

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

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


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

 
  TH1F *h0 = (TH1F*)betaSpec->Get("C9B0");
  TH1F *h1 = (TH1F*)betaSpec->Get("C9B1");
  TH1F *h2 = (TH1F*)betaSpec->Get("C9B2");
  TH1F *h3 = (TH1F*)betaSpec->Get("C9B3");
  TH1F *h4 = (TH1F*)betaSpec->Get("C9B4");
  TH1F *h5 = (TH1F*)betaSpec->Get("C9B5");
  

  for( i=0 ; i<nEvents ; i++ ) {


    if(RandFlat::shoot()<0.0001) {
      cosTheta = RandFlat::shoot(-1, 1);
      azimuth = RandFlat::shoot( 2.0*M_PI );
      theta = acos(cosTheta);
      p0.setRThetaPhi(1, acos(cosTheta), azimuth);

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


    else if(RandFlat::shoot()>=0.0001 && RandFlat::shoot()<=0.0017) {
      cosTheta = RandFlat::shoot(-1, 1);
      azimuth = RandFlat::shoot( 2.0*M_PI );
      theta = acos(cosTheta);
      p0.setRThetaPhi(1, acos(cosTheta), azimuth);

      double ke = h4->GetRandom();
      double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
      p0=p0*pe*keV;
      fprintf(stream, "1\n");
      fprintf(stream, "1\t-11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
    
    }
    else if(RandFlat::shoot()>=0.0017 && RandFlat::shoot()<0.0607) { //beta & gamma 
      cosTheta = RandFlat::shoot(-1, 1);
      azimuth = RandFlat::shoot( 2.0*M_PI );
      theta = acos(cosTheta);
      p0.setRThetaPhi(1, acos(cosTheta), azimuth);

      double ke = h3->GetRandom();
      double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
             p0=p0*pe*keV;

      fprintf(stream, "1\n");
      fprintf(stream, "1\t-11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
    }
    else if(RandFlat::shoot()>=0.0607 && RandFlat::shoot()<0.1187) {                  // beta & gamma
      cosTheta = RandFlat::shoot(-1, 1);
      azimuth = RandFlat::shoot( 2.0*M_PI );
      theta = acos(cosTheta);
      p0.setRThetaPhi(1, acos(cosTheta), azimuth);
      
      double ke = h2->GetRandom();
      double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);

      p0=p0*pe*keV;
      fprintf(stream, "1\n");
      fprintf(stream, "1\t-11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
    }

  else if(RandFlat::shoot()>=0.1187 && RandFlat::shoot()<0.4227) {                  // beta & gamma
      cosTheta = RandFlat::shoot(-1, 1);
      azimuth = RandFlat::shoot( 2.0*M_PI );
      theta = acos(cosTheta);
      p0.setRThetaPhi(1, acos(cosTheta), azimuth);

      double ke = h1->GetRandom();
      double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);

      p0=p0*pe*keV;
      fprintf(stream, "1\n");
      fprintf(stream, "1\t-11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
    }

   else {                  // beta & gamma
      cosTheta = RandFlat::shoot(-1, 1);
      azimuth = RandFlat::shoot( 2.0*M_PI );
      theta = acos(cosTheta);
      p0.setRThetaPhi(1, acos(cosTheta), azimuth);
  
      double ke = h0->GetRandom();
      double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);

      p0=p0*pe*keV;
      fprintf(stream, "1\n");
      fprintf(stream, "1\t-11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.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("C9.exe [-seed seed] [-o outputfilename] [-n nevents]\n");

}

---