In This Package:

AmCNeutronGen.cc

Go to the documentation of this file.

// This is a HepEVT generator of AmC neutron 
// // Modified from K-40 generator by Raymond Tsang, h0237236@hkusua.hku.hk, 30 Aug, 2006.

//\author Jianglai Liu
//jianglai.liu@sjtu.edu.cn, 16 Feb., 2011
#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 "TH1.h"
#include "TH2.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
  Hep3Vector p1; // the gamma momentum vector. Unit GeV/c
  Hep3Vector p2; // 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, "# Ransom seed for generator = %ld.\n", rseed);
  

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

  //cout<<specfilename<<endl;

  TFile  *rootFile = TFile::Open(specfilename.c_str());
  
  if(!rootFile)
    cout << "Can not open "<<specfilename.c_str()<<". "
         << "Make sure SITEROOT or AmC_neutron.root is set to it's location."
         << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/AmC/" << endl;
  if(!(rootFile->IsOpen()))
    cout << "Can not open  AmC_neutron.root. "
         << "Make sure SITEROOT is set to it's location."
         << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/AmC/" << endl;
  
  TH2F *h2 = (TH2F*)rootFile->Get("hEnergyAngle");
  
  double Tn(0);
  double cosTheta(0), theta(0), azimuth(0); // angles used
  Hep3Vector pn; // neutron momentum vector. Unit GeV/c
  for( int i=0 ; i<nEvents ; i++ ) {
    h2->GetRandom2(cosTheta, Tn);
    Tn *= MeV;
    azimuth = RandFlat::shoot( 2.0*M_PI );
    double momentum = sqrt(Tn*Tn + 2*Tn*neutron_mass_c2);
    
    pn.setRThetaPhi(momentum, acos(cosTheta), azimuth);
    
    fprintf(stream, "1\n");
    fprintf(stream, "1\t2112 0 0 %e %e %e %e\n", pn.x()/GeV, pn.y()/GeV, pn.z()/GeV, neutron_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("AmCNeutron.exe [-seed seed] [-o outputfilename] [-n nevents]\n");
        
}

---