Jpp  18.3.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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

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 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1514
Q
Q(UTCMax_s-UTCMin_s)-livetime_s
JPHYSICS::EErad_t
static const JRadiationSource_t EErad_t
Definition: JRadiation.hh:510
E
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
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::DENSITY_SEA_WATER
static const double DENSITY_SEA_WATER
Fixed environment values.
Definition: JPhysics/JConstants.hh:44
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JTOOLS::n
const int n
Definition: JPolint.hh:786
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JMATH::pow
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
JPHYSICS::Brems_t
static const JRadiationSource_t Brems_t
Definition: JRadiation.hh:511
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
N
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
Z
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
Definition: JAlignDetector.sh:81
JTOOLS::j
int j
Definition: JPolint.hh:792
JPHYSICS::GNrad_t
static const JRadiationSource_t GNrad_t
Definition: JRadiation.hh:512
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
debug
int debug
debug level
Definition: archive-put-wiki-detectors.sh:92
RANGE
then fatal Invalid detector identifier $DETECTOR_ID fi set_variable RUNSETUPID typeset a RANGE RANGE[1]
Definition: JCalibratePMT.sh:62
JPHYSICS::JGeaneWater::getA
virtual double getA() const override
Get energy loss constant.
Definition: JGeane.hh:231
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Generated by   doxygen 1.8.5