Jpp  debug
the software that should make you happy
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:69
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1714
JPHYSICS::JGeaneWater::getA
virtual double getA() const override
Get energy loss constant.
Definition: JGeane.hh:231
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
JPHYSICS::GNrad_t
static const JRadiationSource_t GNrad_t
Definition: JRadiation.hh:514
JPHYSICS::DENSITY_SEA_WATER
static const double DENSITY_SEA_WATER
Fixed environment values.
Definition: JPhysics/JConstants.hh:44
JPHYSICS::Brems_t
static const JRadiationSource_t Brems_t
Definition: JRadiation.hh:513
JPHYSICS::H
static const double H
Planck constant [eV s].
Definition: JPhysics/JConstants.hh:29
JPHYSICS::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
JPHYSICS::EErad_t
static const JRadiationSource_t EErad_t
Definition: JRadiation.hh:512
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:786
JTOOLS::j
int j
Definition: JPolint.hh:792
std
Definition: JSTDTypes.hh:14
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.9.1