| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

src

C9DecayGen.cc File Reference

#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 C9DecayGen.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)

Function Documentation

void ProcessArgs ( int  argc,
char **  argv,
long *  rseed,
char **  outfilename,
unsigned int *  nevents 
)

Definition at line 192 of file C9DecayGen.cc.

00193                                                                             {
00194         int i;
00195         for( i=1 ; i<argc ; i++ ) {
00196                 if( strcmp(argv[i], "-seed")==0 ) {
00197                         i++;
00198                         sscanf(argv[i], "%ld", rseed);
00199                                                   }
00200                  else if( strcmp(argv[i], "-o")==0 ) {
00201                         i++;
00202                         *outfilename = new char[strlen(argv[i]) +1];
00203                         strcpy(*outfilename, argv[i]);
00204                                                      } 
00205                 else if( strcmp(argv[i], "-n")==0 ) { 
00206                         i++;
00207                         sscanf(argv[i], "%ud", nevents);
00208                                                     }  
00209                 else {
00210                         Usage();
00211                         exit(0);
00212                      }
00213                                    }
00214 }

void Usage (  ) 

Definition at line 216 of file C9DecayGen.cc.

00216              {
00217         printf("C9.exe [-seed seed] [-o outputfilename] [-n nevents]\n");
00218 
00219 }

int main ( int  argc,
char **  argv 
)

Definition at line 30 of file C9DecayGen.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 
00037 // some more detailed processes will be added later when necessary. 
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/C9/
00057     std::string specfilename = "C9BetaSpectrum.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/C9/";
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 C9BetaSpectrum.root. "
00075                << "Make sure SITEROOT or C9BetaSpectrum.root is set to it's location." 
00076                << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/C9/" << endl;
00077     if(!(betaSpec->IsOpen()))
00078         cout << "Can not open C9BetaSpectrum.root "
00079                << "Make sure SITEROOT is set to it's location."
00080                << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/C9/" << endl;
00081 
00082 
00083   double cosTheta, theta, azimuth; // angles used
00084   unsigned int i;
00085 
00086  
00087   TH1F *h0 = (TH1F*)betaSpec->Get("C9B0");
00088   TH1F *h1 = (TH1F*)betaSpec->Get("C9B1");
00089   TH1F *h2 = (TH1F*)betaSpec->Get("C9B2");
00090   TH1F *h3 = (TH1F*)betaSpec->Get("C9B3");
00091   TH1F *h4 = (TH1F*)betaSpec->Get("C9B4");
00092   TH1F *h5 = (TH1F*)betaSpec->Get("C9B5");
00093   
00094 
00095   for( i=0 ; i<nEvents ; i++ ) {
00096 
00097 
00098     if(RandFlat::shoot()<0.0001) {
00099       cosTheta = RandFlat::shoot(-1, 1);
00100       azimuth = RandFlat::shoot( 2.0*M_PI );
00101       theta = acos(cosTheta);
00102       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00103 
00104       double ke = h5->GetRandom();
00105       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00106       p0=p0*pe*keV;
00107       fprintf(stream, "1\n");
00108       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 );
00109     
00110     }
00111 
00112 
00113     else if(RandFlat::shoot()>=0.0001 && RandFlat::shoot()<=0.0017) {
00114       cosTheta = RandFlat::shoot(-1, 1);
00115       azimuth = RandFlat::shoot( 2.0*M_PI );
00116       theta = acos(cosTheta);
00117       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00118 
00119       double ke = h4->GetRandom();
00120       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00121       p0=p0*pe*keV;
00122       fprintf(stream, "1\n");
00123       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 );
00124     
00125     }
00126     else if(RandFlat::shoot()>=0.0017 && RandFlat::shoot()<0.0607) { //beta & gamma 
00127       cosTheta = RandFlat::shoot(-1, 1);
00128       azimuth = RandFlat::shoot( 2.0*M_PI );
00129       theta = acos(cosTheta);
00130       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00131 
00132       double ke = h3->GetRandom();
00133       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00134              p0=p0*pe*keV;
00135 
00136       fprintf(stream, "1\n");
00137       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 );
00138     }
00139     else if(RandFlat::shoot()>=0.0607 && RandFlat::shoot()<0.1187) {                  // beta & gamma
00140       cosTheta = RandFlat::shoot(-1, 1);
00141       azimuth = RandFlat::shoot( 2.0*M_PI );
00142       theta = acos(cosTheta);
00143       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00144       
00145       double ke = h2->GetRandom();
00146       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00147 
00148       p0=p0*pe*keV;
00149       fprintf(stream, "1\n");
00150       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 );
00151     }
00152 
00153   else if(RandFlat::shoot()>=0.1187 && RandFlat::shoot()<0.4227) {                  // beta & gamma
00154       cosTheta = RandFlat::shoot(-1, 1);
00155       azimuth = RandFlat::shoot( 2.0*M_PI );
00156       theta = acos(cosTheta);
00157       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00158 
00159       double ke = h1->GetRandom();
00160       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00161 
00162       p0=p0*pe*keV;
00163       fprintf(stream, "1\n");
00164       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 );
00165     }
00166 
00167    else {                  // beta & gamma
00168       cosTheta = RandFlat::shoot(-1, 1);
00169       azimuth = RandFlat::shoot( 2.0*M_PI );
00170       theta = acos(cosTheta);
00171       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00172   
00173       double ke = h0->GetRandom();
00174       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00175 
00176       p0=p0*pe*keV;
00177       fprintf(stream, "1\n");
00178       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 );
00179     } 
00180 
00181  
00182   }
00183 
00184 
00185   return 0;
00186 }

| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |
Generated on Mon Apr 11 20:06:43 2011 for C9 by doxygen 1.4.7