Jpp test-rotations-old-533-g2bdbdb559
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 <memory>
#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 "JPhysics/JSeaWater.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 55 of file JRadiation.cc.

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