33int main(
int argc,
char* argv[])
44 JParser<> zap(
"Example program to histogram shower energy.");
52 catch(
const exception &error) {
53 FATAL(error.what() << endl);
59 TH2* h0 =
new TH2D(
"total", NULL, 24, 2.0, 10.0, 20, -5.0, 0.0);
60 TH2* h2 = (TH2*) h0->Clone(
"h2");
62 NOTICE(
"Setting up radiation tables... " << flush);
69 const JRadiation hydrogen ( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
70 const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
71 const JRadiation chlorine (17.0, 35.0, 40, 0.01, 0.1, 0.1);
72 const JRadiation sodium (11.0, 23.0, 40, 0.01, 0.1, 0.1);
79 radiation.push_back(tuple(
new JRadiationSource(11, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), EErad_t), (TH2D*) h0->Clone(
"[eerad O]" )));
80 radiation.push_back(tuple(
new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t), (TH2D*) h0->Clone(
"[eerad Cl]")));
81 radiation.push_back(tuple(
new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), EErad_t), (TH2D*) h0->Clone(
"[eerad H]" )));
82 radiation.push_back(tuple(
new JRadiationSource(14, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t), (TH2D*) h0->Clone(
"[eerad Na]" )));
84 radiation.push_back(tuple(
new JRadiationSource(21, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Brems_t), (TH2D*) h0->Clone(
"[Brems O]" )));
85 radiation.push_back(tuple(
new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t), (TH2D*) h0->Clone(
"[Brems Cl]")));
86 radiation.push_back(tuple(
new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Brems_t), (TH2D*) h0->Clone(
"[Brems H]" )));
87 radiation.push_back(tuple(
new JRadiationSource(24, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t), (TH2D*) h0->Clone(
"[Brems Na]" )));
89 radiation.push_back(tuple(
new JRadiationSource(31, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), GNrad_t), (TH2D*) h0->Clone(
"[gnrad O]" )));
90 radiation.push_back(tuple(
new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t), (TH2D*) h0->Clone(
"[gnrad Cl]")));
91 radiation.push_back(tuple(
new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), GNrad_t), (TH2D*) h0->Clone(
"[gnrad H]" )));
92 radiation.push_back(tuple(
new JRadiationSource(34, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t), (TH2D*) h0->Clone(
"[gnrad Na]" )));
96 for (
int i = 1; i <= h0->GetXaxis()->GetNbins(); ++i) {
98 const double x = h0->GetXaxis()->GetBinCenter(i);
99 const double E = pow(10.0, x);
104 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
108 const double li = p->first->getInverseInteractionLength(E);
109 const double Es = p->first->getEnergyOfShower(E);
110 const double y = log10(Es/E);
112 p->second->Fill(x, y, W);
114 h2->Fill(x, y, li * W);