Jpp 19.3.0-rc.1
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/JQuantile.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
60 string outputFile;
61 int numberOfPoints;
62 string option;
63 int debug;
64
65 try {
66
67 JParser<> zap("Example program to histogram radiation cross sections, shower energy, range and b(E).");
68
69 zap['o'] = make_field(outputFile) = "radiation.root";
70 zap['n'] = make_field(numberOfPoints) = 0;
71 zap['O'] = make_field(option) = LENGTH, ENERGY, RANGE, ELOSS;
72 zap['d'] = make_field(debug) = 0;
73
74 zap(argc, argv);
75 }
76 catch(const exception &error) {
77 FATAL(error.what() << endl);
78 }
79
80
81 using namespace JPP;
82
83
84 TFile out(outputFile.c_str(), "recreate");
85
86 TH1* h0 = NULL; // master histogram
87
88 if (option == LENGTH) {
89 h0 = new TH1D("total", NULL, 160, 2.0, 10.0); // interaction length [m]
90 }
91
92 if (option == ENERGY) {
93 h0 = new TH1D("total", NULL, 24, 2.0, 10.0); // <E> [GeV]
94 }
95
96 NOTICE("Setting up radiation tables... " << flush);
97
98 typedef pair<JSharedPointer<JRadiationInterface>, TH1*> tuple;
99 typedef vector<tuple> ntuple;
100
101 ntuple radiation;
102
103 const JRadiation hydrogen ( 1.0, 1.0, 40, 0.01, 0.1, 0.1);
104 const JRadiation oxygen ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
105 const JRadiation chlorine (17.0, 35.0, 40, 0.01, 0.1, 0.1);
106 const JRadiation sodium (11.0, 23.0, 40, 0.01, 0.1, 0.1);
107
108 JSharedPointer<JRadiation> Hydrogen (new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
109 JSharedPointer<JRadiation> Oxygen (new JRadiationFunction(oxygen, 300, 0.2, 1.0e11));
110 JSharedPointer<JRadiation> Chlorine (new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
111 JSharedPointer<JRadiation> Sodium (new JRadiationFunction(sodium, 300, 0.2, 1.0e11));
112
113 radiation.push_back(tuple(new JRadiationSource(11, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), EErad_t), clone(h0, "[eerad O]" )));
114 radiation.push_back(tuple(new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t), clone(h0, "[eerad Cl]")));
115 radiation.push_back(tuple(new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), EErad_t), clone(h0, "[eerad H]" )));
116 radiation.push_back(tuple(new JRadiationSource(14, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t), clone(h0, "[eerad Na]" )));
117
118 radiation.push_back(tuple(new JRadiationSource(21, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), Brems_t), clone(h0, "[Brems O]" )));
119 radiation.push_back(tuple(new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t), clone(h0, "[Brems Cl]")));
120 radiation.push_back(tuple(new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), Brems_t), clone(h0, "[Brems H]" )));
121 radiation.push_back(tuple(new JRadiationSource(24, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t), clone(h0, "[Brems Na]" )));
122
123 radiation.push_back(tuple(new JRadiationSource(31, Oxygen, DENSITY_SEA_WATER * JSeaWater::O(), GNrad_t), clone(h0, "[gnrad O]" )));
124 radiation.push_back(tuple(new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t), clone(h0, "[gnrad Cl]")));
125 radiation.push_back(tuple(new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(), GNrad_t), clone(h0, "[gnrad H]" )));
126 radiation.push_back(tuple(new JRadiationSource(34, Sodium, DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t), clone(h0, "[gnrad Na]" )));
127
128 NOTICE("OK" << endl);
129
130
131 if (option == LENGTH) {
132
133 for (int i = 1; i <= h0->GetNbinsX(); ++i) {
134
135 const double x = h0->GetBinCenter(i);
136 const double E = pow(10.0, x);
137
138 STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);
139
140 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
141
142 const double li = p->first->getInverseInteractionLength(E);
143
144 p->second->Fill(x, 1.0/li);
145 }
146 }
147 STATUS(endl);
148 }
149
150
151 if (option == ENERGY) {
152
153 for (int i = 1; i <= h0->GetNbinsX(); ++i) {
154
155 const double x = h0->GetBinCenter(i);
156 const double E = pow(10.0, x);
157
158 STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);
159
160 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
161
162 JQuantile Q;
163
164 for (int n = 0; n != numberOfPoints; ++n) {
165 Q.put(p->first->getEnergyOfShower(E));
166 }
167
168 p->second->SetBinContent(i, Q.getMean());
169 //p->second->SetBinError (i, Q.getSTDev());
170 }
171 }
172 STATUS(endl);
173 }
174
175
176 if (option == RANGE) {
177
178 TH1D* Ra = new TH1D("R[analytical]", NULL, 12, 2.0, 8.0); // Range [m]
179 TH1D* Rb = new TH1D("R[numerical]", NULL, 12, 2.0, 8.0); // Range [m]
180
181 for (int i = 1; i <= Ra->GetNbinsX(); ++i) {
182
183 const double x = Ra->GetBinCenter(i);
184 const double E0 = pow(10.0, x);
185 const double Z = gWater(E0);
186
187 Ra->SetBinContent(i, Z*1e-3);
188 }
189
190 for (int j = 1; j <= Rb->GetNbinsX(); ++j) {
191
192 const double x = Rb->GetBinCenter(j);
193 const double E0 = pow(10.0, x);
194
195 STATUS("Energy: " << SCIENTIFIC(8,1) << E0 << '\r'); DEBUG(endl);
196
197 JQuantile Q;
198
199 for (int n = 0; n != numberOfPoints; ++n) {
200
201 double E = E0;
202 double z = 0.0;
203
204 while (E >= 0.5) {
205
206 const int N = radiation.size();
207
208 double li[N]; // inverse interaction lengths
209 double ls = 1.0e-5; // minimal total inverse interaction length (100 km)^-1
210
211 for (int i = 0; i != N; ++i) {
212 ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
213 }
214
215 double dz = min(gRandom->Exp(1.0) / ls, gWater(E));
216
217 E -= gWater.getA() * dz;
218
219 double Es = 0.0; // shower energy [GeV]
220 double y = gRandom->Uniform(ls);
221
222 for (int i = 0; i != N; ++i) {
223
224 y -= li[i];
225
226 if (y < 0.0) {
227 Es = radiation[i].first->getEnergyOfShower(E); // shower energy [GeV]
228 break;
229 }
230 }
231
232 z += dz;
233 E -= Es;
234 }
235
236 Q.put(z);
237 }
238
239 Rb->SetBinContent(j, Q.getMean() * 1e-3);
240 //Rb->SetBinError (j, Q.getSTDev() * 1e-3);
241 }
242 STATUS(endl);
243 }
244
245
246 if (option == ELOSS) {
247
248 TH1D* hb = new TH1D("hb", NULL, 12, 2.0, 8.0); // b(E)
249
250 for (int j = 1; j <= hb->GetNbinsX(); ++j) {
251
252 const double x = hb->GetBinCenter(j);
253 const double E = pow(10.0, x);
254
255 STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);
256
257 JQuantile Q;
258
259 const int N = radiation.size();
260
261 double li[N]; // inverse interaction lengths
262 double ls = 1.0e-5; // minimal total inverse interaction length (100 km)^-1
263
264 for (int i = 0; i != N; ++i) {
265 ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
266 }
267
268 for (int i = 0; i != numberOfPoints; ++i) {
269
270 double Es = 0.0; // shower energy [GeV]
271 double y = gRandom->Uniform(ls);
272
273 for (int i = 0; i != N; ++i) {
274
275 y -= li[i];
276
277 if (y < 0.0) {
278 Es = radiation[i].first->getEnergyOfShower(E); // shower energy [GeV]
279 break;
280 }
281 }
282
283 Q.put(1e3*ls*Es/E);
284 }
285
286 hb->SetBinContent(j, Q.getMean());
287 //hb->SetBinError (j, Q.getSTDev());
288 }
289 STATUS(endl);
290 }
291
292 delete h0;
293
294 out.Write();
295 out.Close();
296}
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:43
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::JQuantile
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
JTOOLS::JQuantile::put
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133
JTOOLS::JQuantile::getMean
double getMean() const
Get mean value.
Definition JQuantile.hh:252
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488
Generated by doxygen 1.12.0