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

In This Package:

src

Be11DecayGen.cc

Go to the documentation of this file.
00001 // This is a HepEVT generator of Be11 decay (including betas & gammas). 
00002 // Modified from K-40 generator by Raymond Tsang, h0237236@hkusua.hku.hk, 30 Aug, 2006.
00003 //
00004 // //\author Fengpeng An
00005 // //anfp@ihep.ac.cn, 8 Nov., 2009
00006 #include <cstdlib>
00007 #include <cstdio>
00008 #include <cstring>
00009 #include <iostream>
00010 #include <math.h>
00011 
00012 #include <CLHEP/Vector/ThreeVector.h>
00013 #include <CLHEP/Random/Randomize.h>
00014 #include <CLHEP/Units/PhysicalConstants.h>
00015 
00016 #include "TFile.h"
00017 #include "TH1F.h"
00018 #include "TCanvas.h"
00019 #include "TRandom.h"
00020 #include "TMath.h"
00021 
00022 using namespace std;
00023 using namespace CLHEP;
00024 
00025 void ProcessArgs(int argc, char** argv, long* rseed, char** outfilename, 
00026                  unsigned int* nevents );
00027 void Usage();
00028 
00029 
00030 int main(int argc, char** argv) {
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   Hep3Vector p3; // the gamma momentum vector. Unit GeV/c
00039   Hep3Vector p4; // the gamma momentum vector. Unit GeV/c
00040   Hep3Vector p5; // the gamma momentum vector. Unit GeV/c
00041   Hep3Vector p6; // the gamma momentum vector. Unit GeV/c
00042   
00043   FILE* stream = stdout;
00044   if( outFilename!=NULL ) {
00045     stream = fopen(outFilename, "w");
00046     if( stream==NULL ) {
00047       printf("Please enter a valid filename.\n");
00048       Usage();
00049       exit(0);
00050     }
00051   }
00052   HepRandom::setTheSeed(rseed);
00053    gRandom->SetSeed(rseed);
00054       
00055   //start by printing some information to comment lines
00056   fprintf(stream, "# File generated by %s.\n", argv[0]);
00057   fprintf(stream, "# Random seed for generator = %ld.\n", rseed);
00058   
00059 
00060 // The root file used by this generator should be located at 
00061 // $SITEROOT/dybgaudi/Generators/RadioActivity/Be11/
00062     std::string specfilename = "Be11BetaSpectrum.root";
00063     const char* prefix = getenv("SITEROOT");
00064     if(prefix) {
00065         if(strlen(prefix)>0) {
00066             std::string newname = prefix;
00067             newname += "/dybgaudi/Generators/RadioActivity/Be11/";
00068             newname += specfilename;
00069             specfilename = newname;
00070         }
00071     }
00072 
00073 
00074 //cout<<specfilename<<endl;
00075 
00076   TFile  *betaSpec = TFile::Open(specfilename.c_str());
00077 
00078     if(!betaSpec)
00079         cout << "Can not open Be11BetaSpectrum.root. "
00080                << "Make sure SITEROOT or Be11BetaSpectrum.root is set to it's location."
00081                << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/Be11/" << endl;
00082     if(!(betaSpec->IsOpen()))
00083         cout << "Can not open Be11BetaSpectrum.root "
00084                << "Make sure SITEROOT is set to it's location."
00085                << "rootfile location is at $SITEROOT/dybgaudi/Generators/RadioActivity/Be11/" << endl;
00086 
00087 
00088   double cosTheta, theta, azimuth; // angles used
00089   unsigned int i;
00090 
00091   TH1F *h0 = (TH1F*)betaSpec->Get("Be11B0");
00092   TH1F *h1 = (TH1F*)betaSpec->Get("Be11B1");
00093   TH1F *h2 = (TH1F*)betaSpec->Get("Be11B2");
00094   TH1F *h3 = (TH1F*)betaSpec->Get("Be11B3");
00095   TH1F *h4 = (TH1F*)betaSpec->Get("Be11B4");
00096   TH1F *h5 = (TH1F*)betaSpec->Get("Be11B5");
00097   TH1F *h6 = (TH1F*)betaSpec->Get("Be11B6");
00098   
00099   for( i=0 ; i<nEvents ; i++ ) {
00100 
00101     if(RandFlat::shoot()<=0.031) { //beta & gamma 
00102       cosTheta = RandFlat::shoot(-1, 1);
00103       azimuth = RandFlat::shoot( 2.0*M_PI );
00104       theta = acos(cosTheta);
00105       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00106 
00107       double ke = h6->GetRandom();
00108       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00109       p0=p0*pe*keV;
00110       fprintf(stream, "1\n");
00111       fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00112     
00113     }                
00114 
00115 
00116     else if(RandFlat::shoot()>0.031 && RandFlat::shoot()<= 0.071) { //beta & gamma 
00117       cosTheta = RandFlat::shoot(-1, 1);
00118       azimuth = RandFlat::shoot( 2.0*M_PI );
00119       theta = acos(cosTheta);
00120       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00121 
00122       cosTheta = RandFlat::shoot(-1, 1);
00123       azimuth = RandFlat::shoot( 2.0*M_PI );
00124       theta = acos(cosTheta);
00125       p1.setRThetaPhi(1, acos(cosTheta), azimuth);
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       cosTheta = RandFlat::shoot(-1, 1);
00133       azimuth = RandFlat::shoot( 2.0*M_PI );
00134       theta = acos(cosTheta);
00135       p3.setRThetaPhi(1, acos(cosTheta), azimuth);
00136 
00137       cosTheta = RandFlat::shoot(-1, 1);
00138       azimuth = RandFlat::shoot( 2.0*M_PI );
00139       theta = acos(cosTheta);
00140       p4.setRThetaPhi(1, acos(cosTheta), azimuth);
00141 
00142       cosTheta = RandFlat::shoot(-1, 1);
00143       azimuth = RandFlat::shoot( 2.0*M_PI );
00144       theta = acos(cosTheta);
00145       p5.setRThetaPhi(1, acos(cosTheta), azimuth);
00146 
00147       double ke = h5->GetRandom();
00148       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00149              p0=p0*pe*keV;
00150 
00151            if(RandFlat::shoot()<=0.462)
00152               {
00153                  p1*=7.97473*MeV;
00154 
00155                  fprintf(stream, "2\n");
00156                  fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00157                  fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV );
00158               }    
00159            else if(RandFlat::shoot()>0.462 && RandFlat::shoot()<= 0.994)
00160               { 
00161                 p2*=5.85147*MeV;
00162                 p3*=4.4439*MeV;
00163 
00164                 fprintf(stream, "3\n");
00165                 fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00166                 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p2.x()/GeV, p2.y()/GeV, p2.z()/GeV );
00167                 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p3.x()/GeV, p3.y()/GeV, p3.z()/GeV );
00168                }
00169              
00170           else
00171              { 
00172                 p4*=0.692*MeV;
00173                 p5*=7.2829*MeV;
00174 
00175                if(RandFlat::shoot()<=0.870)
00176                   {
00177                   fprintf(stream, "3\n");
00178                   fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00179                   fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p4.x()/GeV, p4.y()/GeV, p4.z()/GeV );
00180                   fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p5.x()/GeV, p5.y()/GeV, p5.z()/GeV );
00181                   }
00182 
00183                else
00184                  {
00185                  fprintf(stream, "2\n");
00186                  fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00187                  fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p4.x()/GeV, p4.y()/GeV, p4.z()/GeV );
00188                  }
00189 
00190              }
00191 
00192     }
00193     else if(RandFlat::shoot()>0.071 && RandFlat::shoot()<=0.136) {                  // beta & gamma
00194       
00195       cosTheta = RandFlat::shoot(-1, 1);
00196       azimuth = RandFlat::shoot( 2.0*M_PI );
00197       theta = acos(cosTheta);
00198       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00199 
00200       cosTheta = RandFlat::shoot(-1, 1);
00201       azimuth = RandFlat::shoot( 2.0*M_PI );
00202       theta = acos(cosTheta);
00203       p1.setRThetaPhi(1, acos(cosTheta), azimuth);
00204 
00205       cosTheta = RandFlat::shoot(-1, 1);
00206       azimuth = RandFlat::shoot( 2.0*M_PI );
00207       theta = acos(cosTheta);
00208       p2.setRThetaPhi(1, acos(cosTheta), azimuth);
00209       
00210       cosTheta = RandFlat::shoot(-1, 1);
00211       azimuth = RandFlat::shoot( 2.0*M_PI );
00212       theta = acos(cosTheta);
00213       p3.setRThetaPhi(1, acos(cosTheta), azimuth);
00214 
00215       cosTheta = RandFlat::shoot(-1, 1);
00216       azimuth = RandFlat::shoot( 2.0*M_PI );
00217       theta = acos(cosTheta);
00218       p4.setRThetaPhi(1, acos(cosTheta), azimuth);
00219 
00220       cosTheta = RandFlat::shoot(-1, 1);
00221       azimuth = RandFlat::shoot( 2.0*M_PI );
00222       theta = acos(cosTheta);
00223       p5.setRThetaPhi(1, acos(cosTheta), azimuth);
00224 
00225       cosTheta = RandFlat::shoot(-1, 1);
00226       azimuth = RandFlat::shoot( 2.0*M_PI );
00227       theta = acos(cosTheta);
00228       p6.setRThetaPhi(1, acos(cosTheta), azimuth);
00229 
00230       double ke = h4->GetRandom();
00231       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00232       p0=p0*pe*keV;
00233            
00234            if(RandFlat::shoot()<=0.675)
00235               {
00236                  p1*=6.78981*MeV;
00237 
00238                   fprintf(stream, "2\n");
00239                  fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00240                  fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV );
00241               }    
00242            else if(RandFlat::shoot()>0.675 && RandFlat::shoot()<= 0.96)
00243               { 
00244                 p2*=4.6659*MeV;
00245                 p3*=4.4439*MeV;
00246 
00247                 fprintf(stream, "3\n");
00248                 fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00249                 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p2.x()/GeV, p2.y()/GeV, p2.z()/GeV );
00250                 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p3.x()/GeV, p3.y()/GeV, p3.z()/GeV );
00251               }
00252 
00253            else
00254              {
00255                if(RandFlat::shoot()<=0.856)
00256                {
00257                  p4*=1.77131*MeV;
00258                  p5*=5.01898*MeV;
00259                  fprintf(stream, "3\n");
00260                  fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00261                  fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p4.x()/GeV, p4.y()/GeV, p4.z()/GeV );
00262                  fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p5.x()/GeV, p5.y()/GeV, p5.z()/GeV );
00263                
00264                }
00265                else 
00266                {
00267                  p4*=1.77131*MeV;
00268                  p6*=2.89530*MeV;
00269                  fprintf(stream, "3\n");
00270                  fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00271                  fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p4.x()/GeV, p4.y()/GeV, p4.z()/GeV );
00272                  fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p6.x()/GeV, p6.y()/GeV, p6.z()/GeV );
00273 
00274                }
00275            
00276              }
00277 
00278 
00279 
00280     }
00281     
00282     else if(RandFlat::shoot()>0.136 && RandFlat::shoot()<=0.1388)
00283     {
00284 
00285       cosTheta = RandFlat::shoot(-1, 1);
00286       azimuth = RandFlat::shoot( 2.0*M_PI );
00287       theta = acos(cosTheta);
00288       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00289 
00290       cosTheta = RandFlat::shoot(-1, 1);
00291       azimuth = RandFlat::shoot( 2.0*M_PI );
00292       theta = acos(cosTheta);
00293       p1.setRThetaPhi(1, acos(cosTheta), azimuth);
00294 
00295       cosTheta = RandFlat::shoot(-1, 1);
00296       azimuth = RandFlat::shoot( 2.0*M_PI );
00297       theta = acos(cosTheta);
00298       p2.setRThetaPhi(1, acos(cosTheta), azimuth);
00299 
00300       cosTheta = RandFlat::shoot(-1, 1);
00301       azimuth = RandFlat::shoot( 2.0*M_PI );
00302       theta = acos(cosTheta);
00303       p3.setRThetaPhi(1, acos(cosTheta), azimuth);
00304 
00305       double ke = h3->GetRandom();
00306       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00307       p0=p0*pe*keV;
00308 
00309            if(RandFlat::shoot()<=0.856)
00310               {
00311                  p1*=5.01898*MeV;
00312 
00313                  fprintf(stream, "2\n");
00314                  fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00315                  fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV );
00316               }    
00317            else 
00318               { 
00319                 p2*=2.8953*MeV;
00320                 p3*=4.4439*MeV;
00321 
00322                 fprintf(stream, "3\n");
00323                 fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00324                 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p2.x()/GeV, p2.y()/GeV, p2.z()/GeV );
00325                 fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p3.x()/GeV, p3.y()/GeV, p3.z()/GeV );
00326                }
00327     
00328     }
00329 
00330     else if(RandFlat::shoot()>0.1388 && RandFlat::shoot()<=0.13934)
00331     {
00332 
00333       cosTheta = RandFlat::shoot(-1, 1);
00334       azimuth = RandFlat::shoot( 2.0*M_PI );
00335       theta = acos(cosTheta);
00336       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00337 
00338       cosTheta = RandFlat::shoot(-1, 1);
00339       azimuth = RandFlat::shoot( 2.0*M_PI );
00340       theta = acos(cosTheta);
00341       p1.setRThetaPhi(1, acos(cosTheta), azimuth);
00342 
00343       double ke = h2->GetRandom();
00344       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00345       p0=p0*pe*keV;
00346       p1*=4.4439*MeV;
00347       fprintf(stream, "2\n");
00348       fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00349       fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV );
00350     }
00351 
00352     else if(RandFlat::shoot()>0.13934 && RandFlat::shoot()<=0.453339)
00353     {
00354 
00355       cosTheta = RandFlat::shoot(-1, 1);
00356       azimuth = RandFlat::shoot( 2.0*M_PI );
00357       theta = acos(cosTheta);
00358       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00359 
00360       cosTheta = RandFlat::shoot(-1, 1);
00361       azimuth = RandFlat::shoot( 2.0*M_PI );
00362       theta = acos(cosTheta);
00363       p1.setRThetaPhi(1, acos(cosTheta), azimuth);
00364 
00365       double ke = h1->GetRandom();
00366       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00367       p0=p0*pe*keV;
00368       p1*=2.12447*MeV;
00369       fprintf(stream, "2\n");
00370       fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00371       fprintf(stream, "1\t22 0 0 %e %e %e 0.0\n", p1.x()/GeV, p1.y()/GeV, p1.z()/GeV );
00372      }
00373 
00374     else 
00375     {
00376       cosTheta = RandFlat::shoot(-1, 1);
00377       azimuth = RandFlat::shoot( 2.0*M_PI );
00378       theta = acos(cosTheta);
00379       p0.setRThetaPhi(1, acos(cosTheta), azimuth);
00380 
00381       double ke = h0->GetRandom();
00382       double pe = sqrt(ke*ke + 2*ke*electron_mass_c2/keV);
00383       p0=p0*pe*keV;
00384       fprintf(stream, "1\n");
00385       fprintf(stream, "1\t11 0 0 %e %e %e %e\n", p0.x()/GeV, p0.y()/GeV, p0.z()/GeV, electron_mass_c2/GeV );
00386     }
00387   }
00388   return 0;
00389 }
00390 
00391 
00392 
00393 
00394 
00395 void ProcessArgs(int argc, char** argv, long *rseed, 
00396                                  char** outfilename, unsigned int* nevents) {
00397         int i;
00398         for( i=1 ; i<argc ; i++ ) {
00399                 if( strcmp(argv[i], "-seed")==0 ) {
00400                         i++;
00401                         sscanf(argv[i], "%ld", rseed);
00402                                                   }
00403                  else if( strcmp(argv[i], "-o")==0 ) {
00404                         i++;
00405                         *outfilename = new char[strlen(argv[i]) +1];
00406                         strcpy(*outfilename, argv[i]);
00407                                                      } 
00408                 else if( strcmp(argv[i], "-n")==0 ) { 
00409                         i++;
00410                         sscanf(argv[i], "%ud", nevents);
00411                                                     }  
00412                 else {
00413                         Usage();
00414                         exit(0);
00415                      }
00416                                    }
00417 }
00418 
00419 void Usage() {
00420         printf("Be11.exe [-seed seed] [-o outputfilename] [-n nevents]\n");
00421 
00422 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |
Generated on Mon Apr 11 20:06:52 2011 for Be11 by doxygen 1.4.7