00001 #include <iostream>
00002 #include <stdlib.h>
00003 #include <utility>
00004
00005
00006 #include <CLHEP/Random/JamesRandom.h>
00007
00008
00009 #include "../CrElectron.h"
00010
00011 typedef double G4double;
00012
00013 const int BIN_NUM = 1024;
00014 const int LOOP = int(1e+6);
00015
00016 const G4double energy_min = 0.1;
00017 const G4double energy_max = 100.0;
00018
00019
00020 int main(int argc, char** argv)
00021 {
00022 using namespace std;
00023 const double GeV=1, rad=1;
00024 int hist[BIN_NUM];
00025 for (int i = 0; i < BIN_NUM; i++) hist[i] = 0;
00026
00027 CrElectron src;
00028
00029 HepRandomEngine* engine = new HepJamesRandom;
00030
00031 int loop = (argc > 1) ? (int)atof(*(argv + argc - 1)) : LOOP;
00032 {for (int i = 0; i < loop; i++){
00033 src.selectComponent(engine);
00034 G4double energy = src.energySrc(engine) / GeV;
00035 std::pair<double,double> dir = src.dir(energy, engine);
00036
00037 G4double theta = acos(dir.first);
00038 G4double phi = dir.second * rad;
00039 hist[int(BIN_NUM * (energy - energy_min) / (energy_max - energy_min))]++;
00040 if (i % 10000 == 1) cerr << '\015' << i << ": " << energy << "...";
00041 }}
00042 cerr << "\n";
00043
00044 {for (int i = 0; i < BIN_NUM; i++){
00045 cout << ((i + 0.5)/BIN_NUM) * (energy_max - energy_min) + energy_min << " "
00046 << hist[i] << "\n";
00047 }}
00048
00049 delete engine;
00050
00051 return 0;
00052 }