Jpp  19.1.0
the software that should make you happy
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 
11 #include "JPhysics/JRadiation.hh"
14 #include "JPhysics/JGeane.hh"
15 #include "JPhysics/JConstants.hh"
16 
17 #include "JSirene/JSeaWater.hh"
18 #include "JLang/JSharedPointer.hh"
19 
20 #include "Jeep/JPrint.hh"
21 #include "Jeep/JParser.hh"
22 #include "Jeep/JMessage.hh"
23 
24 
25 namespace {
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  */
56 int main(int argc, char* argv[])
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 }
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:69
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
Utility class to parse command line options.
Definition: JParser.hh:1698
virtual double getA() const override
Get energy loss constant.
Definition: JGeane.hh:231
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
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
static const JRadiationSource_t GNrad_t
Definition: JRadiation.hh:514
static const double DENSITY_SEA_WATER
Fixed environment values.
static const JRadiationSource_t Brems_t
Definition: JRadiation.hh:513
static const double H
Planck constant [eV s].
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
static const JRadiationSource_t EErad_t
Definition: JRadiation.hh:512
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
const int n
Definition: JPolint.hh:786
int j
Definition: JPolint.hh:792
Definition: JSTDTypes.hh:14
Auxiliary data structure to list files in directory.
Definition: JFilesystem.hh:20
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