#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"Include dependency graph for N12DecayGen.cc:
Go to the source code of this file.
Namespaces | |
| namespace | std |
| namespace | CLHEP |
Functions | |
| void | ProcessArgs (int argc, char **argv, long *rseed, char **outfilename, unsigned int *nevents) |
| void | Usage () |
| int | main (int argc, char **argv) |
| void ProcessArgs | ( | int | argc, | |
| char ** | argv, | |||
| long * | rseed, | |||
| char ** | outfilename, | |||
| unsigned int * | nevents | |||
| ) |
Definition at line 166 of file N12DecayGen.cc.
00167 { 00168 int i; 00169 for( i=1 ; i<argc ; i++ ) { 00170 if( strcmp(argv[i], "-seed")==0 ) { 00171 i++; 00172 sscanf(argv[i], "%ld", rseed); 00173 } 00174 else if( strcmp(argv[i], "-o")==0 ) { 00175 i++; 00176 *outfilename = new char[strlen(argv[i]) +1]; 00177 strcpy(*outfilename, argv[i]); 00178 } 00179 else if( strcmp(argv[i], "-n")==0 ) { 00180 i++; 00181 sscanf(argv[i], "%ud", nevents); 00182 } 00183 else { 00184 Usage(); 00185 exit(0); 00186 } 00187 } 00188 }
| void Usage | ( | ) |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 31 of file N12DecayGen.cc.
00031 { 00032 long rseed = 0; 00033 char* outFilename = NULL; 00034 unsigned int nEvents = 1000000000; // a billion. Default to something big for piping 00035 ProcessArgs(argc, argv, &rseed, &outFilename, &nEvents); 00036 Hep3Vector p0; // the beta momentum vector. Unit GeV/c 00037 Hep3Vector p1; // the gamma momentum vector. Unit GeV/c 00038 Hep3Vector p2; // the gamma momentum vector. Unit GeV/c 00039 00040 FILE* stream = stdout; 00041 if( outFilename!=NULL ) { 00042 stream = fopen(outFilename, "w"); 00043 if( stream==NULL ) { 00044 printf("Please enter a valid filename.\n"); 00045 Usage(); 00046 exit(0); 00047 } 00048 } 00049 HepRandom::setTheSeed(rseed); 00050 gRandom->SetSeed(rseed); 00051 //start by printing some information to comment lines 00052 fprintf(stream, "# File generated by %s.\n", argv[0]); 00053 fprintf(stream, "# Random seed for generator = %ld.\n", rseed); 00054 00055 00056 // The root file used by this generator should be located at 00057 // $SITEROOT/dybgaudi/Generators/RadioActivity/N12/ 00058 std::string specfilename = "N12BetaSpectrum.root"; 00059 const char* prefix = getenv("SITEROOT"); 00060 if(prefix) { 00061 if(strlen(prefix)>0) { 00062 std::string newname = prefix; 00063 newname += "/dybgaudi/Generators/RadioActivity/N12/"; 00064 newname += specfilename; 00065 specfilename = newname; 00066 } 00067 } 00068 00069 00070 //cout<<specfilename<<endl; 00071 00072 TFile *betaSpec = TFile::Open(specfilename.c_str()); 00073 00074 if(!betaSpec) 00075 cout << "Can not open N12BetaSpectrum.root. " 00076 << "Make sure SITEROOT or N12BetaSpectrum.root is set to it's location." 00077 << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/N12/" << endl; 00078 if(!(betaSpec->IsOpen())) 00079 cout << "Can not open N12BetaSpectrum.root " 00080 << "Make sure SITEROOT is set to it's location." 00081 << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/N12/" << endl; 00082 00083 double cosTheta, theta, azimuth; // angles used 00084 unsigned int i; 00085 00086 TH1F *h0 = (TH1F*)betaSpec->Get("N12B0"); 00087 TH1F *h1 = (TH1F*)betaSpec->Get("N12B1"); 00088 TH1F *h2 = (TH1F*)betaSpec->Get("N12B2"); 00089 00090 for( i=0 ; i<nEvents ; i++ ) { 00091 00092 if(RandFlat::shoot()<0.027) { //beta & gamma 00093 cosTheta = RandFlat::shoot(-1, 1); 00094 azimuth = RandFlat::shoot( 2.0*M_PI ); 00095 theta = acos(cosTheta); 00096 p0.setRThetaPhi(1, acos(cosTheta), azimuth); 00097 00098 cosTheta = RandFlat::shoot(-1, 1); 00099 azimuth = RandFlat::shoot( 2.0*M_PI ); 00100 theta = acos(cosTheta); 00101 p1.setRThetaPhi(1, acos(cosTheta), azimuth); 00102 00103 cosTheta = RandFlat::shoot(-1, 1); 00104 azimuth = RandFlat::shoot( 2.0*M_PI ); 00105 theta = acos(cosTheta); 00106 p2.setRThetaPhi(1, acos(cosTheta), azimuth); 00107 00108 double ke = h2->GetRandom(); 00109 double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV); 00110 00111 p0=p0*pe*keV; 00112 p1*=3.21483*MeV; 00113 p2*=4.43803*MeV; 00114 fprintf(stream, "3\n"); 00115 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 ); 00116 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV ); 00117 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p2.x()/GeV, p2.y()/GeV, p2.z()/GeV ); 00118 00119 } 00120 else if(RandFlat::shoot()>=0.027 && RandFlat::shoot()<0.046) { // beta & gamma 00121 cosTheta = RandFlat::shoot(-1, 1); 00122 azimuth = RandFlat::shoot( 2.0*M_PI ); 00123 theta = acos(cosTheta); 00124 p0.setRThetaPhi(1, acos(cosTheta), azimuth); 00125 00126 00127 cosTheta = RandFlat::shoot(-1, 1); 00128 azimuth = RandFlat::shoot( 2.0*M_PI ); 00129 theta = acos(cosTheta); 00130 p2.setRThetaPhi(1, acos(cosTheta), azimuth); 00131 00132 00133 double ke = h1->GetRandom(); 00134 double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV); 00135 00136 p0=p0*pe*keV; 00137 p2*=4.43803*MeV; 00138 fprintf(stream, "2\n"); 00139 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 ); 00140 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p2.x()/GeV, p2.y()/GeV, p2.z()/GeV ); 00141 } 00142 00143 else { 00144 cosTheta = RandFlat::shoot(-1, 1); 00145 azimuth = RandFlat::shoot( 2.0*M_PI ); 00146 theta = acos(cosTheta); 00147 p0.setRThetaPhi(1, acos(cosTheta), azimuth); 00148 00149 double ke = h0->GetRandom(); 00150 double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV); 00151 p0=p0*pe*keV; 00152 fprintf(stream, "1\n"); 00153 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 ); 00154 00155 } 00156 } 00157 00158 00159 return 0; 00160 }
1.4.7