In This Package:

C10DecayGen.cc

Go to the documentation of this file.

// This is a HepEVT generator of C10 decay (resulting positrons & gammas). 
// // 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 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/C10/
    std::string specfilename = "C10BetaSpectrum.root";
    const char* prefix = getenv("SITEROOT");
    if(prefix) {
        if(strlen(prefix)>0) {
            std::string newname = prefix;
            newname += "/dybgaudi/Generators/RadioActivity/C10/";
            newname += specfilename;
            specfilename = newname;
        }
    }


//cout<<specfilename<<endl;

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

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


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

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

    if(RandFlat::shoot()<0.01465) {  

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

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

      cosTheta = RandFlat::shoot(-1, 1);
      azimuth = RandFlat::shoot( 2.0*M_PI );
      theta = acos(cosTheta);
      p2.setRThetaPhi(1, acos(cosTheta), azimuth);
      
      double ke = h2->GetRandom();
      double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);

      p0=p0*pe*keV;
      p1*=1.0217*MeV;
      p2*=0.7183*MeV;
      fprintf(stream, "3\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 );
      fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV );
      fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p2.x()/GeV, p2.y()/GeV, p2.z()/GeV );

    }

    else {

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

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

      p0=p0*pe*keV;
      p1*=0.7183*MeV;
      fprintf(stream, "2\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 );
      fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/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("C10.exe [-seed seed] [-o outputfilename] [-n nevents]\n");

}

---