#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 C10DecayGen.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 151 of file C10DecayGen.cc.
00152 { 00153 int i; 00154 for( i=1 ; i<argc ; i++ ) { 00155 if( strcmp(argv[i], "-seed")==0 ) { 00156 i++; 00157 sscanf(argv[i], "%ld", rseed); 00158 } 00159 else if( strcmp(argv[i], "-o")==0 ) { 00160 i++; 00161 *outfilename = new char[strlen(argv[i]) +1]; 00162 strcpy(*outfilename, argv[i]); 00163 } 00164 else if( strcmp(argv[i], "-n")==0 ) { 00165 i++; 00166 sscanf(argv[i], "%ud", nevents); 00167 } 00168 else { 00169 Usage(); 00170 exit(0); 00171 } 00172 } 00173 }
| void Usage | ( | ) |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 30 of file C10DecayGen.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, "# Ransom seed for generator = %ld.\n", rseed); 00053 00054 00055 // The root file used by this generator should be located at 00056 // $SITEROOT/dybgaudi/Generators/RadioActivity/C10/ 00057 std::string specfilename = "C10BetaSpectrum.root"; 00058 const char* prefix = getenv("SITEROOT"); 00059 if(prefix) { 00060 if(strlen(prefix)>0) { 00061 std::string newname = prefix; 00062 newname += "/dybgaudi/Generators/RadioActivity/C10/"; 00063 newname += specfilename; 00064 specfilename = newname; 00065 } 00066 } 00067 00068 00069 //cout<<specfilename<<endl; 00070 00071 TFile *betaSpec = TFile::Open(specfilename.c_str()); 00072 00073 if(!betaSpec) 00074 cout << "Can not open C10BetaSpectrum.root. " 00075 << "Make sure SITEROOT or C10BetaSpectrum.root is set to it's location." 00076 << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/C10/" << endl; 00077 if(!(betaSpec->IsOpen())) 00078 cout << "Can not open C10BetaSpectrum.root " 00079 << "Make sure SITEROOT is set to it's location." 00080 << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/C10/" << endl; 00081 00082 00083 double cosTheta, theta, azimuth; // angles used 00084 unsigned int i; 00085 00086 TH1F *h1 = (TH1F*)betaSpec->Get("C10B1"); 00087 TH1F *h2 = (TH1F*)betaSpec->Get("C10B2"); 00088 00089 for( i=0 ; i<nEvents ; i++ ) { 00090 00091 if(RandFlat::shoot()<0.01465) { 00092 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*=1.0217*MeV; 00113 p2*=0.7183*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 00121 else { 00122 00123 cosTheta = RandFlat::shoot(-1, 1); 00124 azimuth = RandFlat::shoot( 2.0*M_PI ); 00125 theta = acos(cosTheta); 00126 p0.setRThetaPhi(1, acos(cosTheta), azimuth); 00127 00128 cosTheta = RandFlat::shoot(-1, 1); 00129 azimuth = RandFlat::shoot( 2.0*M_PI ); 00130 theta = acos(cosTheta); 00131 p1.setRThetaPhi(1, acos(cosTheta), azimuth); 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 p1*=0.7183*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", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV ); 00141 } 00142 } 00143 00144 return 0; 00145 }
1.4.7