Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JRadiation.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4
5#include "TROOT.h"
6#include "TFile.h"
7#include "TH1D.h"
8
9#include "JTools/JQuantile.hh"
10
14#include "JPhysics/JGeane.hh"
16
17#include "JSirene/JSeaWater.hh"
19
20#include "Jeep/JPrint.hh"
21#include "Jeep/JParser.hh"
22#include "Jeep/JMessage.hh"
23
24
25namespace {
26
27 static const std::string LENGTH = "length"; //!< calculation of interaction length
28 static const std::string ENERGY = "energy"; //!< calculation of shower energy
29 static const std::string RANGE = "range"; //!< calculation of muon range
30 static const std::string ELOSS = "eloss"; //!< calculation of b parameter in dE/dx formula
31
32 /**
33 * Clone master histogram.
34 *
35 * \param h1 pointer to master histogram
36 * \param buffer name of clone
37 * \return pointer to cloned histogram
38 */
39 inline TH1* clone(TH1* h1, const char* buffer)
40 {
41 if (h1 != NULL)
42 return (TH1*) h1->Clone(buffer);
43 else
44 return NULL;
45 }
46}
47
48
49/**
50 * \file
51 *
52 * Example program to histogram radiation cross sections, shower energy, range and b(E)
53 * as a function of the muon energy.
54 * \author mdejong
55 */
56int main(int argc, char* argv[])
57{
58 using namespace std;
59 using namespace JPP;
60
61 string outputFile;
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
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 JQuantile 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 JQuantile 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 JQuantile 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}
string outputFile
Energy loss of muon.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
I/O formatting auxiliaries.
int main(int argc, char *argv[])
Definition JRadiation.cc:56
Muon radiative cross sections.
int numberOfPoints
Definition JResultPDF.cc:22
The template JSharedPointer class can be used to share a pointer to an object.
Utility class to parse command line options.
Definition JParser.hh:1698
Fast implementation of class JRadiation.
Implementation for calculation of inverse interaction length and shower energy.
Auxiliary class for the calculation of the muon radiative cross sections.
Definition JRadiation.hh:36
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary data structure to list files in directory.
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133
double getMean() const
Get mean value.
Definition JQuantile.hh:252
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488