27 static const std::string LENGTH =
"length";
28 static const std::string ENERGY =
"energy";
29 static const std::string RANGE =
"range";
30 static const std::string ELOSS =
"eloss";
39 inline TH1* clone(TH1*
h1,
const char* buffer)
42 return (TH1*) h1->Clone(buffer);
56 int main(
int argc,
char* argv[])
67 JParser<> zap(
"Example program to histogram radiation cross sections, shower energy, range and b(E).");
71 zap[
'O'] =
make_field(option) = LENGTH, ENERGY, RANGE, ELOSS;
76 catch(
const exception &error) {
77 FATAL(error.what() << endl);
88 if (option == LENGTH) {
89 h0 =
new TH1D(
"total", NULL, 160, 2.0, 10.0);
92 if (option == ENERGY) {
93 h0 =
new TH1D(
"total", NULL, 24, 2.0, 10.0);
96 NOTICE(
"Setting up radiation tables... " << flush);
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);
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));
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]")));
116 radiation.push_back(tuple(
new JRadiationSource(14, Sodium,
DENSITY_SEA_WATER * JSeaWater::Na(),
EErad_t), clone(h0,
"[eerad Na]" )));
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]")));
121 radiation.push_back(tuple(
new JRadiationSource(24, Sodium,
DENSITY_SEA_WATER * JSeaWater::Na(),
Brems_t), clone(h0,
"[Brems Na]" )));
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]")));
126 radiation.push_back(tuple(
new JRadiationSource(34, Sodium,
DENSITY_SEA_WATER * JSeaWater::Na(),
GNrad_t), clone(h0,
"[gnrad Na]" )));
131 if (option == LENGTH) {
133 for (
int i = 1; i <= h0->GetNbinsX(); ++i) {
135 const double x = h0->GetBinCenter(i);
136 const double E =
pow(10.0, x);
140 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
142 const double li = p->first->getInverseInteractionLength(E);
144 p->second->Fill(x, 1.0/li);
151 if (option == ENERGY) {
153 for (
int i = 1; i <= h0->GetNbinsX(); ++i) {
155 const double x = h0->GetBinCenter(i);
156 const double E =
pow(10.0, x);
160 for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
165 Q.put(p->first->getEnergyOfShower(E));
168 p->second->SetBinContent(i, Q.getMean());
176 if (option == RANGE) {
178 TH1D* Ra =
new TH1D(
"R[analytical]", NULL, 12, 2.0, 8.0);
179 TH1D* Rb =
new TH1D(
"R[numerical]", NULL, 12, 2.0, 8.0);
181 for (
int i = 1; i <= Ra->GetNbinsX(); ++i) {
183 const double x = Ra->GetBinCenter(i);
184 const double E0 =
pow(10.0, x);
187 Ra->SetBinContent(i, Z*1e-3);
190 for (
int j = 1;
j <= Rb->GetNbinsX(); ++
j) {
192 const double x = Rb->GetBinCenter(
j);
193 const double E0 =
pow(10.0, x);
206 const int N = radiation.size();
211 for (
int i = 0; i !=
N; ++i) {
212 ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
215 double dz = min(gRandom->Exp(1.0) / ls,
gWater(E));
220 double y = gRandom->Uniform(ls);
222 for (
int i = 0; i !=
N; ++i) {
227 Es = radiation[i].first->getEnergyOfShower(E);
239 Rb->SetBinContent(
j, Q.getMean() * 1e-3);
246 if (option == ELOSS) {
248 TH1D* hb =
new TH1D(
"hb", NULL, 12, 2.0, 8.0);
250 for (
int j = 1;
j <= hb->GetNbinsX(); ++
j) {
252 const double x = hb->GetBinCenter(
j);
253 const double E =
pow(10.0, x);
259 const int N = radiation.size();
264 for (
int i = 0; i !=
N; ++i) {
265 ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
271 double y = gRandom->Uniform(ls);
273 for (
int i = 0; i !=
N; ++i) {
278 Es = radiation[i].first->getEnergyOfShower(E);
286 hb->SetBinContent(
j, Q.getMean());
Utility class to parse command line options.
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
static const JRadiationSource_t EErad_t
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
static const double H
Planck constant [eV s].
then for HISTOGRAM in h0 h1
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
static const double DENSITY_SEA_WATER
Fixed environment values.
Muon radiative cross sections.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
T pow(const T &x, const double y)
Power .
static const JRadiationSource_t Brems_t
General purpose messaging.
Utility class to parse command line options.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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
static const JRadiationSource_t GNrad_t
Auxiliary data structure for floating point format specification.
virtual double getA() const override
Get energy loss constant.