Jpp 20.0.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JRadiation.cc File Reference

Example program to histogram radiation cross sections, shower energy, range and b(E) as a function of the muon energy. More...

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "JTools/JStats.hh"
#include "JPhysics/JRadiation.hh"
#include "JPhysics/JRadiationFunction.hh"
#include "JPhysics/JRadiationSource.hh"
#include "JPhysics/JGeane.hh"
#include "JPhysics/JConstants.hh"
#include "JSirene/JSeaWater.hh"
#include "JLang/JSharedPointer.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Detailed Description

Example program to histogram radiation cross sections, shower energy, range and b(E) as a function of the muon energy.

Author
mdejong

Definition in file JRadiation.cc.

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 56 of file JRadiation.cc.

57{
58 using namespace std;
59 using namespace JPP;
60
61 string outputFile;
62 int numberOfPoints;
63 string option;
64 int debug;
65
66 try {
67
68 JParser<> zap("Example program to histogram radiation cross sections, shower energy, range and b(E).");
69
70 zap['o'] = make_field(outputFile) = "radiation.root";
71 zap['n'] = make_field(numberOfPoints) = 0;
72 zap['O'] = make_field(option) = LENGTH, ENERGY, RANGE, ELOSS;
73 zap['d'] = make_field(debug) = 0;
74
75 zap(argc, argv);
76 }
77 catch(const exception &error) {
78 FATAL(error.what() << endl);
79 }
80
81
82 TFile out(outputFile.c_str(), "recreate");
83
84 TH1* h0 = NULL; // master histogram
85
86 if (option == LENGTH) {
87 h0 = new TH1D("total", NULL, 160, 2.0, 10.0); // interaction length [m]
88 }
89
90 if (option == ENERGY) {
91 h0 = new TH1D("total", NULL, 24, 2.0, 10.0); // <E> [GeV]
92 }
93
94 NOTICE("Setting up radiation tables... " << flush);
95
96 typedef pair<JSharedPointer<JRadiationInterface>, TH1*> tuple;
97 typedef vector<tuple> ntuple;
98
99 ntuple radiation;
100
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);
105
106 JSharedPointer<JRadiation> Hydrogen (new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
107 JSharedPointer<JRadiation> Oxygen (new JRadiationFunction(oxygen, 300, 0.2, 1.0e11));
108 JSharedPointer<JRadiation> Chlorine (new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
109 JSharedPointer<JRadiation> Sodium (new JRadiationFunction(sodium, 300, 0.2, 1.0e11));
110
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]" )));
115
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]" )));
120
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]" )));
125
126 NOTICE("OK" << endl);
127
128
129 if (option == LENGTH) {
130
131 for (int i = 1; i <= h0->GetNbinsX(); ++i) {
132
133 const double x = h0->GetBinCenter(i);
134 const double E = pow(10.0, x);
135
136 STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);
137
138 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
139
140 const double li = p->first->getInverseInteractionLength(E);
141
142 p->second->Fill(x, 1.0/li);
143 }
144 }
145 STATUS(endl);
146 }
147
148
149 if (option == ENERGY) {
150
151 for (int i = 1; i <= h0->GetNbinsX(); ++i) {
152
153 const double x = h0->GetBinCenter(i);
154 const double E = pow(10.0, x);
155
156 STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);
157
158 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
159
160 JStats Q;
161
162 for (int n = 0; n != numberOfPoints; ++n) {
163 Q.put(p->first->getEnergyOfShower(E));
164 }
165
166 p->second->SetBinContent(i, Q.getMean());
167 //p->second->SetBinError (i, Q.getSTDev());
168 }
169 }
170 STATUS(endl);
171 }
172
173
174 if (option == RANGE) {
175
176 TH1D* Ra = new TH1D("R[analytical]", NULL, 12, 2.0, 8.0); // Range [m]
177 TH1D* Rb = new TH1D("R[numerical]", NULL, 12, 2.0, 8.0); // Range [m]
178
179 for (int i = 1; i <= Ra->GetNbinsX(); ++i) {
180
181 const double x = Ra->GetBinCenter(i);
182 const double E0 = pow(10.0, x);
183 const double Z = gWater(E0);
184
185 Ra->SetBinContent(i, Z*1e-3);
186 }
187
188 for (int j = 1; j <= Rb->GetNbinsX(); ++j) {
189
190 const double x = Rb->GetBinCenter(j);
191 const double E0 = pow(10.0, x);
192
193 STATUS("Energy: " << SCIENTIFIC(8,1) << E0 << '\r'); DEBUG(endl);
194
195 JStats Q;
196
197 for (int n = 0; n != numberOfPoints; ++n) {
198
199 double E = E0;
200 double z = 0.0;
201
202 while (E >= 0.5) {
203
204 const int N = radiation.size();
205
206 double li[N]; // inverse interaction lengths
207 double ls = 1.0e-5; // minimal total inverse interaction length (100 km)^-1
208
209 for (int i = 0; i != N; ++i) {
210 ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
211 }
212
213 double dz = min(gRandom->Exp(1.0) / ls, gWater(E));
214
215 E -= gWater.getA() * dz;
216
217 double Es = 0.0; // shower energy [GeV]
218 double y = gRandom->Uniform(ls);
219
220 for (int i = 0; i != N; ++i) {
221
222 y -= li[i];
223
224 if (y < 0.0) {
225 Es = radiation[i].first->getEnergyOfShower(E); // shower energy [GeV]
226 break;
227 }
228 }
229
230 z += dz;
231 E -= Es;
232 }
233
234 Q.put(z);
235 }
236
237 Rb->SetBinContent(j, Q.getMean() * 1e-3);
238 //Rb->SetBinError (j, Q.getSTDev() * 1e-3);
239 }
240 STATUS(endl);
241 }
242
243
244 if (option == ELOSS) {
245
246 TH1D* hb = new TH1D("hb", NULL, 12, 2.0, 8.0); // b(E)
247
248 for (int j = 1; j <= hb->GetNbinsX(); ++j) {
249
250 const double x = hb->GetBinCenter(j);
251 const double E = pow(10.0, x);
252
253 STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);
254
255 JStats Q;
256
257 const int N = radiation.size();
258
259 double li[N]; // inverse interaction lengths
260 double ls = 1.0e-5; // minimal total inverse interaction length (100 km)^-1
261
262 for (int i = 0; i != N; ++i) {
263 ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
264 }
265
266 for (int i = 0; i != numberOfPoints; ++i) {
267
268 double Es = 0.0; // shower energy [GeV]
269 double y = gRandom->Uniform(ls);
270
271 for (int i = 0; i != N; ++i) {
272
273 y -= li[i];
274
275 if (y < 0.0) {
276 Es = radiation[i].first->getEnergyOfShower(E); // shower energy [GeV]
277 break;
278 }
279 }
280
281 Q.put(1e3*ls*Es/E);
282 }
283
284 hb->SetBinContent(j, Q.getMean());
285 //hb->SetBinError (j, Q.getSTDev());
286 }
287 STATUS(endl);
288 }
289
290 delete h0;
291
292 out.Write();
293 out.Close();
294}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
numberOfPoints
int numberOfPoints
Definition JResultPDF.cc:22
JLANG::JSharedPointer
The template JSharedPointer class can be used to share a pointer to an object.
Definition JSharedPointer.hh:31
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JPHYSICS::JRadiationFunction
Fast implementation of class JRadiation.
Definition JRadiationFunction.hh:34
JPHYSICS::JRadiationSource
Implementation for calculation of inverse interaction length and shower energy.
Definition JRadiationSource.hh:93
JPHYSICS::JRadiation
Auxiliary class for the calculation of the muon radiative cross sections.
Definition JRadiation.hh:36
JMATH::pow
T pow(const T &x, const double y)
Power .
Definition JMath.hh:97
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JTOOLS::n
const int n
Definition JPolint.hh:791
JTOOLS::j
int j
Definition JPolint.hh:801
JSYSTEM::ls
Auxiliary data structure to list files in directory.
Definition JFilesystem.hh:20
JTOOLS::JStats
Auxiliary data structure for running average and standard deviation.
Definition JStats.hh:44
JTOOLS::JStats::put
void put(const double x, const double w=1.0)
Put value.
Definition JStats.hh:119
JTOOLS::JStats::getMean
double getMean() const
Get mean value.
Definition JStats.hh:234
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488
Generated by doxygen 1.12.0