#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 B12DecayGen.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 B12DecayGen.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 30 of file B12DecayGen.cc.
00030 { 00031 long rseed = 0; 00032 char* outFilename = NULL; 00033 unsigned int nEvents = 1000000000; // a billion. Default to something big for piping 00034 ProcessArgs(argc, argv, &rseed, &outFilename, &nEvents); 00035 Hep3Vector p0; // the beta momentum vector. Unit GeV/c 00036 Hep3Vector p1; // the gamma momentum vector. Unit GeV/c 00037 Hep3Vector p2; // the gamma momentum vector. Unit GeV/c 00038 00039 FILE* stream = stdout; 00040 if( outFilename!=NULL ) { 00041 stream = fopen(outFilename, "w"); 00042 if( stream==NULL ) { 00043 printf("Please enter a valid filename.\n"); 00044 Usage(); 00045 exit(0); 00046 } 00047 } 00048 HepRandom::setTheSeed(rseed); 00049 gRandom->SetSeed(rseed); 00050 //start by printing some information to comment lines 00051 fprintf(stream, "# File generated by %s.\n", argv[0]); 00052 fprintf(stream, "# Random seed for generator = %ld.\n", rseed); 00053 00054 // The root file used by this generator should be located at 00055 // $SITEROOT/dybgaudi/Generators/RadioActivity/B12/ 00056 std::string specfilename = "B12BetaSpectrum.root"; 00057 const char* prefix = getenv("SITEROOT"); 00058 if(prefix) { 00059 if(strlen(prefix)>0) { 00060 std::string newname = prefix; 00061 newname += "/dybgaudi/Generators/RadioActivity/B12/"; 00062 newname += specfilename; 00063 specfilename = newname; 00064 } 00065 } 00066 00067 00068 cout<<specfilename<<endl; 00069 00070 TFile *betaSpec = TFile::Open(specfilename.c_str()); 00071 00072 if(!betaSpec) 00073 cout << "Can not open B12BetaSpectrum.root. " 00074 << "Make sure SITEROOT or B12BetaSpectrum.root is set to it's location." 00075 << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/B12/" << endl; 00076 if(!(betaSpec->IsOpen())) 00077 cout << "Can not open B12BetaSpectrum.root " 00078 << "Make sure SITEROOT is set to it's location." 00079 << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/B12/" << endl; 00080 00081 00082 00083 double cosTheta, theta, azimuth; // angles used 00084 unsigned int i; 00085 00086 TH1F *h0 = (TH1F*)betaSpec->Get("B12B0"); 00087 TH1F *h1 = (TH1F*)betaSpec->Get("B12B1"); 00088 TH1F *h2 = (TH1F*)betaSpec->Get("B12B2"); 00089 00090 for( i=0 ; i<nEvents ; i++ ) { 00091 00092 if(RandFlat::shoot()<0.015) { //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\t11 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.015 && RandFlat::shoot()<0.0273) { // 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\t11 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\t11 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