56int main(
int argc,
char* argv[])
68 JParser<> zap(
"Example program to histogram radiation cross sections, shower energy, range and b(E).");
72 zap[
'O'] =
make_field(option) = LENGTH, ENERGY, RANGE, ELOSS;
77 catch(
const exception &error) {
78 FATAL(error.what() << endl);
86 if (option == LENGTH) {
87 h0 =
new TH1D(
"total", NULL, 160, 2.0, 10.0);
90 if (option == ENERGY) {
91 h0 =
new TH1D(
"total", NULL, 24, 2.0, 10.0);
94 NOTICE(
"Setting up radiation tables... " << flush);
101 const JRadiation hydrogen ( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
102 const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
103 const JRadiation chlorine (17.0, 35.0, 40, 0.01, 0.1, 0.1);
104 const JRadiation sodium (11.0, 23.0, 40, 0.01, 0.1, 0.1);
111 radiation.push_back(tuple(
new JRadiationSource(11, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), EErad_t), clone(h0,
"[eerad O]" )));
112 radiation.push_back(tuple(
new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t), clone(h0,
"[eerad Cl]")));
113 radiation.push_back(tuple(
new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), EErad_t), clone(h0,
"[eerad H]" )));
114 radiation.push_back(tuple(
new JRadiationSource(14, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t), clone(h0,
"[eerad Na]" )));
116 radiation.push_back(tuple(
new JRadiationSource(21, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Brems_t), clone(h0,
"[Brems O]" )));
117 radiation.push_back(tuple(
new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t), clone(h0,
"[Brems Cl]")));
118 radiation.push_back(tuple(
new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Brems_t), clone(h0,
"[Brems H]" )));
119 radiation.push_back(tuple(
new JRadiationSource(24, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t), clone(h0,
"[Brems Na]" )));
121 radiation.push_back(tuple(
new JRadiationSource(31, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), GNrad_t), clone(h0,
"[gnrad O]" )));
122 radiation.push_back(tuple(
new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t), clone(h0,
"[gnrad Cl]")));
123 radiation.push_back(tuple(
new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), GNrad_t), clone(h0,
"[gnrad H]" )));
124 radiation.push_back(tuple(
new JRadiationSource(34, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t), clone(h0,
"[gnrad Na]" )));
129 if (option == LENGTH) {
131 for (
int i = 1; i <= h0->GetNbinsX(); ++i) {
133 const double x = h0->GetBinCenter(i);
134 const double E = pow(10.0, x);
138 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
140 const double li = p->first->getInverseInteractionLength(E);
142 p->second->Fill(x, 1.0/li);
149 if (option == ENERGY) {
151 for (
int i = 1; i <= h0->GetNbinsX(); ++i) {
153 const double x = h0->GetBinCenter(i);
154 const double E = pow(10.0, x);
158 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
163 Q.
put(p->first->getEnergyOfShower(E));
166 p->second->SetBinContent(i, Q.
getMean());
174 if (option == RANGE) {
176 TH1D* Ra =
new TH1D(
"R[analytical]", NULL, 12, 2.0, 8.0);
177 TH1D* Rb =
new TH1D(
"R[numerical]", NULL, 12, 2.0, 8.0);
179 for (
int i = 1; i <= Ra->GetNbinsX(); ++i) {
181 const double x = Ra->GetBinCenter(i);
182 const double E0 = pow(10.0, x);
183 const double Z = gWater(E0);
185 Ra->SetBinContent(i, Z*1e-3);
188 for (
int j = 1; j <= Rb->GetNbinsX(); ++j) {
190 const double x = Rb->GetBinCenter(j);
191 const double E0 = pow(10.0, x);
204 const int N = radiation.size();
209 for (
int i = 0; i != N; ++i) {
210 ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
213 double dz = min(gRandom->Exp(1.0) /
ls, gWater(E));
215 E -= gWater.getA() * dz;
218 double y = gRandom->Uniform(
ls);
220 for (
int i = 0; i != N; ++i) {
225 Es = radiation[i].first->getEnergyOfShower(E);
237 Rb->SetBinContent(j, Q.
getMean() * 1e-3);
244 if (option == ELOSS) {
246 TH1D* hb =
new TH1D(
"hb", NULL, 12, 2.0, 8.0);
248 for (
int j = 1; j <= hb->GetNbinsX(); ++j) {
250 const double x = hb->GetBinCenter(j);
251 const double E = pow(10.0, x);
257 const int N = radiation.size();
262 for (
int i = 0; i != N; ++i) {
263 ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
269 double y = gRandom->Uniform(
ls);
271 for (
int i = 0; i != N; ++i) {
276 Es = radiation[i].first->getEnergyOfShower(E);
284 hb->SetBinContent(j, Q.
getMean());