REST-for-Physics  v2.3
Rare Event Searches ToolKit for Physics
TRestGeant4PrimaryGeneratorInfo.cxx
1//
2// Created by lobis on 7/14/2022.
3//
4
5#include "TRestGeant4PrimaryGeneratorInfo.h"
6
7#include <TAxis.h>
8#include <TF1.h>
9#include <TF2.h>
10#include <TMath.h>
11#include <TRestStringOutput.h>
12#include <tinyxml.h>
13
14#include <iostream>
15
16using namespace std;
17using namespace TRestGeant4PrimaryGeneratorTypes;
18
19string TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypesToString(const SpatialGeneratorTypes& type) {
20 switch (type) {
21 case SpatialGeneratorTypes::CUSTOM:
22 return "Custom";
23 case SpatialGeneratorTypes::VOLUME:
24 return "Volume";
25 case SpatialGeneratorTypes::SURFACE:
26 return "Surface";
27 case SpatialGeneratorTypes::POINT:
28 return "Point";
29 case SpatialGeneratorTypes::COSMIC:
30 return "Cosmic";
31 }
32 cout << "TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypesToString - Error - Unknown "
33 "SpatialGeneratorTypes"
34 << endl;
35 exit(1);
36}
37
38SpatialGeneratorTypes TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorTypes(const string& type) {
39 if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::CUSTOM),
40 TString::ECaseCompare::kIgnoreCase)) {
41 return SpatialGeneratorTypes::CUSTOM;
42 } else if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::VOLUME),
43 TString::ECaseCompare::kIgnoreCase)) {
44 return SpatialGeneratorTypes::VOLUME;
45 } else if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::SURFACE),
46 TString::ECaseCompare::kIgnoreCase)) {
47 return SpatialGeneratorTypes::SURFACE;
48 } else if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::POINT),
49 TString::ECaseCompare::kIgnoreCase)) {
50 return SpatialGeneratorTypes::POINT;
51 } else if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::COSMIC),
52 TString::ECaseCompare::kIgnoreCase)) {
53 return SpatialGeneratorTypes::COSMIC;
54 } else {
55 cout << "TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorTypes - Error - Unknown "
56 "SpatialGeneratorTypes: "
57 << type << endl;
58 exit(1);
59 }
60}
61
62string TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorShapesToString(const SpatialGeneratorShapes& type) {
63 switch (type) {
64 case SpatialGeneratorShapes::GDML:
65 return "GDML";
66 case SpatialGeneratorShapes::WALL:
67 return "Wall";
68 case SpatialGeneratorShapes::CIRCLE:
69 return "Circle";
70 case SpatialGeneratorShapes::BOX:
71 return "Box";
72 case SpatialGeneratorShapes::SPHERE:
73 return "Sphere";
74 case SpatialGeneratorShapes::CYLINDER:
75 return "Cylinder";
76 }
77 cout << "TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorShapesToString - Error - Unknown "
78 "SpatialGeneratorShapes"
79 << endl;
80 exit(1);
81}
82
83SpatialGeneratorShapes TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorShapes(const string& type) {
84 if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::GDML),
85 TString::ECaseCompare::kIgnoreCase)) {
86 return SpatialGeneratorShapes::GDML;
87 } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::WALL),
88 TString::ECaseCompare::kIgnoreCase)) {
89 return SpatialGeneratorShapes::WALL;
90 } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::CIRCLE),
91 TString::ECaseCompare::kIgnoreCase)) {
92 return SpatialGeneratorShapes::CIRCLE;
93 } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::BOX),
94 TString::ECaseCompare::kIgnoreCase)) {
95 return SpatialGeneratorShapes::BOX;
96 } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::SPHERE),
97 TString::ECaseCompare::kIgnoreCase)) {
98 return SpatialGeneratorShapes::SPHERE;
99 } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::CYLINDER),
100 TString::ECaseCompare::kIgnoreCase)) {
101 return SpatialGeneratorShapes::CYLINDER;
102
103 } else {
104 cout << "TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorShapes - Error - Unknown "
105 "SpatialGeneratorShapes: "
106 << type << endl;
107 exit(1);
108 }
109}
110
111string TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypesToString(
112 const EnergyDistributionTypes& type) {
113 switch (type) {
114 case EnergyDistributionTypes::TH1D:
115 return "TH1D";
116 case EnergyDistributionTypes::FORMULA:
117 return "Formula";
118 case EnergyDistributionTypes::FORMULA2:
119 return "Formula2";
120 case EnergyDistributionTypes::MONO:
121 return "Mono";
122 case EnergyDistributionTypes::FLAT:
123 return "Flat";
124 case EnergyDistributionTypes::LOG:
125 return "Log";
126 }
127 cout << "TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypesToString - Error - Unknown energy "
128 "distribution type"
129 << endl;
130 exit(1);
131}
132
133EnergyDistributionTypes TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes(
134 const string& type) {
135 if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::TH1D),
136 TString::ECaseCompare::kIgnoreCase)) {
137 return EnergyDistributionTypes::TH1D;
138 } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::FORMULA),
139 TString::ECaseCompare::kIgnoreCase)) {
140 return EnergyDistributionTypes::FORMULA;
141 } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::FORMULA2),
142 TString::ECaseCompare::kIgnoreCase)) {
143 return EnergyDistributionTypes::FORMULA2;
144 } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::MONO),
145 TString::ECaseCompare::kIgnoreCase)) {
146 return EnergyDistributionTypes::MONO;
147 } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::FLAT),
148 TString::ECaseCompare::kIgnoreCase)) {
149 return EnergyDistributionTypes::FLAT;
150 } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::LOG),
151 TString::ECaseCompare::kIgnoreCase)) {
152 return EnergyDistributionTypes::LOG;
153 } else {
154 cout << "TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes - Error - Unknown "
155 "EnergyDistributionTypes: "
156 << type << endl;
157 exit(1);
158 }
159}
160
161string TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToString(
162 const EnergyDistributionFormulas& type) {
163 switch (type) {
164 case EnergyDistributionFormulas::COSMIC_NEUTRONS:
165 return "CosmicNeutrons";
166 case EnergyDistributionFormulas::COSMIC_GAMMAS:
167 return "CosmicGammas";
168 }
169 cout << "TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToString - Error - Unknown energy "
170 "distribution formula"
171 << endl;
172 exit(1);
173}
174
175EnergyDistributionFormulas TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionFormulas(
176 const string& type) {
177 if (TString(type).EqualTo(EnergyDistributionFormulasToString(EnergyDistributionFormulas::COSMIC_NEUTRONS),
178 TString::ECaseCompare::kIgnoreCase)) {
179 return EnergyDistributionFormulas::COSMIC_NEUTRONS;
180 } else if (TString(type).EqualTo(
181 EnergyDistributionFormulasToString(EnergyDistributionFormulas::COSMIC_GAMMAS),
182 TString::ECaseCompare::kIgnoreCase)) {
183 return EnergyDistributionFormulas::COSMIC_GAMMAS;
184 } else {
185 cout << "TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionFormulas - Error - Unknown "
186 "energyDistributionFormulas: "
187 << type << endl;
188 exit(1);
189 }
190}
191
192TF1 TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToRootFormula(
193 const EnergyDistributionFormulas& type) {
194 switch (type) {
195 case EnergyDistributionFormulas::COSMIC_NEUTRONS: {
196 // Formula from https://ieeexplore.ieee.org/document/1369506
197 const char* title = "Cosmic Neutrons at Sea Level";
198 auto distribution =
199 TF1(title,
200 "1.006E-6 * TMath::Exp(-0.3500 * TMath::Power(TMath::Log(x * 1E-3), 2) + 2.1451 * "
201 "TMath::Log(x * 1E-3)) + "
202 "1.011E-3 * TMath::Exp(-0.4106 * TMath::Power(TMath::Log(x * 1E-3), 2) - 0.6670 * "
203 "TMath::Log(x * 1E-3))",
204 1E2, 1E7);
205 distribution.SetNormalized(true); // Normalized, not really necessary
206 distribution.SetTitle(title);
207 distribution.GetXaxis()->SetTitle("Energy (keV)");
208 return distribution;
209 }
210 case EnergyDistributionFormulas::COSMIC_GAMMAS:
211 exit(1);
212 }
213 cout << "TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToRootFormula - Error - Unknown "
214 "energy distribution formula"
215 << endl;
216 exit(1);
217}
218
219string TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypesToString(
220 const AngularDistributionTypes& type) {
221 switch (type) {
222 case AngularDistributionTypes::TH1D:
223 return "TH1D";
224 case AngularDistributionTypes::FORMULA:
225 return "Formula";
226 case AngularDistributionTypes::FORMULA2:
227 return "Formula2";
228 case AngularDistributionTypes::ISOTROPIC:
229 return "Isotropic";
230 case AngularDistributionTypes::FLUX:
231 return "Flux";
232 case AngularDistributionTypes::BACK_TO_BACK:
233 return "Back to back";
234 }
235 cout << "TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypesToString - Error - Unknown angular "
236 "distribution type"
237 << endl;
238 exit(1);
239}
240
241AngularDistributionTypes TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(
242 const string& type) {
243 if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::TH1D),
244 TString::ECaseCompare::kIgnoreCase)) {
245 return AngularDistributionTypes::TH1D;
246 } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::FORMULA),
247 TString::ECaseCompare::kIgnoreCase)) {
248 return AngularDistributionTypes::FORMULA;
249 } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::FORMULA2),
250 TString::ECaseCompare::kIgnoreCase)) {
251 return AngularDistributionTypes::FORMULA2;
252 } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::ISOTROPIC),
253 TString::ECaseCompare::kIgnoreCase)) {
254 return AngularDistributionTypes::ISOTROPIC;
255 } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::FLUX),
256 TString::ECaseCompare::kIgnoreCase)) {
257 return AngularDistributionTypes::FLUX;
258 } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::BACK_TO_BACK),
259 TString::ECaseCompare::kIgnoreCase)) {
260 return AngularDistributionTypes::BACK_TO_BACK;
261 } else {
262 cout << "TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes - Error - Unknown "
263 "AngularDistributionTypes: "
264 << type << endl;
265 exit(1);
266 }
267}
268
269string TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToString(
270 const AngularDistributionFormulas& type) {
271 switch (type) {
272 case AngularDistributionFormulas::COS2:
273 return "Cos2";
274 case AngularDistributionFormulas::COS3:
275 return "Cos3";
276 }
277 cout << "TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToString - Error - Unknown angular "
278 "distribution formula"
279 << endl;
280 exit(1);
281}
282
283AngularDistributionFormulas TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionFormulas(
284 const string& type) {
285 if (TString(type).EqualTo(AngularDistributionFormulasToString(AngularDistributionFormulas::COS2),
286 TString::ECaseCompare::kIgnoreCase)) {
287 return AngularDistributionFormulas::COS2;
288 } else if (TString(type).EqualTo(AngularDistributionFormulasToString(AngularDistributionFormulas::COS3),
289 TString::ECaseCompare::kIgnoreCase)) {
290 return AngularDistributionFormulas::COS3;
291 } else {
292 cout << "TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionFormulas - Error - Unknown "
293 "AngularDistributionFormulas: "
294 << type << endl;
295 exit(1);
296 }
297}
298
299TF1 TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToRootFormula(
300 const AngularDistributionFormulas& type) {
301 switch (type) {
302 case AngularDistributionFormulas::COS2: {
303 auto cos2 = [](double* xs, double* ps) {
304 if (xs[0] >= 0 && xs[0] <= TMath::Pi() / 2) {
305 return TMath::Power(TMath::Cos(xs[0]), 2);
306 }
307 return 0.0;
308 };
309 const char* title = "AngularDistribution: Cos2";
310 auto f = TF1(title, cos2, 0.0, TMath::Pi());
311 f.SetTitle(title);
312 return f;
313 }
314 case AngularDistributionFormulas::COS3: {
315 auto cos3 = [](double* xs, double* ps) {
316 if (xs[0] >= 0 && xs[0] <= TMath::Pi() / 2) {
317 return TMath::Power(TMath::Cos(xs[0]), 3);
318 }
319 return 0.0;
320 };
321 const char* title = "AngularDistribution: Cos3";
322 auto f = TF1(title, cos3, 0.0, TMath::Pi());
323 f.SetTitle(title);
324 return f;
325 }
326 }
327 cout << "TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToRootFormula - Error - Unknown "
328 "angular distribution formula"
329 << endl;
330 exit(1);
331}
332
333string TRestGeant4PrimaryGeneratorTypes::EnergyAndAngularDistributionFormulasToString(
334 const EnergyAndAngularDistributionFormulas& type) {
335 switch (type) {
336 case EnergyAndAngularDistributionFormulas::COSMIC_MUONS:
337 return "CosmicMuons";
338 }
339 cout << "TRestGeant4PrimaryGeneratorTypes::EnergyAndAngularDistributionFormulasToString - Error - "
340 "Unknown energy/angular distribution formula"
341 << endl;
342 exit(1);
343}
344
345EnergyAndAngularDistributionFormulas
346TRestGeant4PrimaryGeneratorTypes::StringToEnergyAndAngularDistributionFormulas(const string& type) {
347 if (TString(type).EqualTo(
348 EnergyAndAngularDistributionFormulasToString(EnergyAndAngularDistributionFormulas::COSMIC_MUONS),
349 TString::ECaseCompare::kIgnoreCase)) {
350 return EnergyAndAngularDistributionFormulas::COSMIC_MUONS;
351 } else {
352 cout << "TRestGeant4PrimaryGeneratorTypes::StringToEnergyAndAngularDistributionFormulas - Error - "
353 "Unknown AngularDistributionFormulas: "
354 << type << endl;
355 exit(1);
356 }
357}
358
359TF2 TRestGeant4PrimaryGeneratorTypes::EnergyAndAngularDistributionFormulasToRootFormula(
360 const EnergyAndAngularDistributionFormulas& type) {
361 switch (type) {
362 case EnergyAndAngularDistributionFormulas::COSMIC_MUONS: {
363 // Guan formula from https://arxiv.org/pdf/1509.06176.pdf
364 // muon rest mass is 105.7 MeV
365 // energy in keV
366 // already integrated in phi (*2pi): formula returns (counts/cm2/s)/keV/rad(theta)
367 const char* title = "Cosmic Muons Energy and Angular";
368 auto f =
369 TF2(title,
370 "2*TMath::Pi()*1E-6*0.14*TMath::Power(x*1E-6*(1.+3.64/"
371 "(x*1E-6*TMath::Power(TMath::Power((TMath::Power(TMath::Cos(y),2)+0.0105212-0.068287*"
372 "TMath::Power(TMath::Cos(y),0.958633)+0.0407253*TMath::Power(TMath::Cos(y),0.817285)"
373 ")/(0.982960),0.5),1.29))),-2.7)*(1./"
374 "(1.+(1.1*x*1E-6*TMath::Power((TMath::Power(TMath::Cos(y),2)+0.0105212-0.068287*TMath::"
375 "Power(TMath::Cos(y),0.958633)+0.0407253*TMath::Power(TMath::Cos(y),0.817285))/"
376 "(0.982960),0.5))/115.)+0.054/"
377 "(1.+(1.1*x*1E-6*TMath::Power((TMath::Power(TMath::Cos(y),2)+0.0105212-0.068287*TMath::"
378 "Power(TMath::Cos(y),0.958633)+0.0407253*TMath::Power(TMath::Cos(y),0.817285))/"
379 "(0.982960),0.5))/850.))*(2.*TMath::Sin(y)*TMath::Pi())",
380 2.0E5, 5.0E9, 0, TMath::Pi() / 2.);
381 f.SetTitle(title);
382 /*
383 * we need to increase the number of bins to get a smooth distribution
384 * this makes the first random generation slow, but it is only a constant factor which does not
385 * matter in long simulations
386 * If we don't include this, the distribution will be very discrete. Max value is 10000
387 */
388 return f;
389 }
390 }
391 cout << "TRestGeant4PrimaryGeneratorTypes::EnergyAndAngularDistributionFormulasToRootFormula - Error - "
392 "Unknown energy and angular distribution formula"
393 << endl;
394 exit(1);
395}
396
397void TRestGeant4PrimaryGeneratorInfo::Print() const {
398 const auto typeEnum = StringToSpatialGeneratorTypes(fSpatialGeneratorType.Data());
399 RESTMetadata << "Generator type: " << SpatialGeneratorTypesToString(typeEnum) << RESTendl;
400
401 if (typeEnum != SpatialGeneratorTypes::POINT) {
402 const auto shapeEnum = StringToSpatialGeneratorShapes(fSpatialGeneratorShape.Data());
403 RESTMetadata << "Generator shape: " << SpatialGeneratorShapesToString(shapeEnum);
404 if (shapeEnum == SpatialGeneratorShapes::GDML) {
405 RESTMetadata << "::" << fSpatialGeneratorFrom << RESTendl;
406 } else {
407 if (shapeEnum == SpatialGeneratorShapes::BOX) {
408 RESTMetadata << ", (length, width, height): ";
409 } else if (shapeEnum == SpatialGeneratorShapes::SPHERE) {
410 RESTMetadata << ", (radius, , ): ";
411 } else if (shapeEnum == SpatialGeneratorShapes::WALL) {
412 RESTMetadata << ", (length, width, ): ";
413 } else if (shapeEnum == SpatialGeneratorShapes::CIRCLE) {
414 RESTMetadata << ", (radius, , ): ";
415 } else if (shapeEnum == SpatialGeneratorShapes::CYLINDER) {
416 RESTMetadata << ", (radius, length, ): ";
417 }
418
419 RESTMetadata << fSpatialGeneratorSize.X() << ", " << fSpatialGeneratorSize.Y() << ", "
420 << fSpatialGeneratorSize.Z() << RESTendl;
421 }
422 }
423 RESTMetadata << "Generator center : (" << fSpatialGeneratorPosition.X() << ","
424 << fSpatialGeneratorPosition.Y() << "," << fSpatialGeneratorPosition.Z() << ") mm"
425 << RESTendl;
426 RESTMetadata << "Generator rotation : (" << fSpatialGeneratorRotationAxis.X() << ","
427 << fSpatialGeneratorRotationAxis.Y() << "," << fSpatialGeneratorRotationAxis.Z()
428 << "), angle: " << fSpatialGeneratorRotationValue * TMath::RadToDeg() << "ยบ" << RESTendl;
429}
TRestGeant4PrimaryGeneratorInfo::fSpatialGeneratorSize
TVector3 fSpatialGeneratorSize
The size of the generator. Stores up to three dimensions according to the shape box: length,...
Definition: TRestGeant4PrimaryGeneratorInfo.h:124
TRestGeant4PrimaryGeneratorInfo::fSpatialGeneratorFrom
TString fSpatialGeneratorFrom
The volume name where the events are generated, in case of volume or surface generator types.
Definition: TRestGeant4PrimaryGeneratorInfo.h:106
TRestGeant4PrimaryGeneratorInfo::fSpatialGeneratorRotationValue
Double_t fSpatialGeneratorRotationValue
degrees of rotation for generator virtual shape along the axis
Definition: TRestGeant4PrimaryGeneratorInfo.h:116
TRestGeant4PrimaryGeneratorInfo::fSpatialGeneratorPosition
TVector3 fSpatialGeneratorPosition
The position of the generator for virtual, and point, generator types.
Definition: TRestGeant4PrimaryGeneratorInfo.h:109
TRestGeant4PrimaryGeneratorInfo::fSpatialGeneratorRotationAxis
TVector3 fSpatialGeneratorRotationAxis
A 3d-std::vector with the angles, measured in degrees, of a XYZ rotation applied to the virtual gener...
Definition: TRestGeant4PrimaryGeneratorInfo.h:113
TRestGeant4PrimaryGeneratorInfo::fSpatialGeneratorShape
TString fSpatialGeneratorShape
Shape of spatial generator (wall, GDML, sphere, etc)
Definition: TRestGeant4PrimaryGeneratorInfo.h:103
TRestGeant4PrimaryGeneratorInfo::fSpatialGeneratorType
TString fSpatialGeneratorType
Type of spatial generator (point, surface, volume, custom)
Definition: TRestGeant4PrimaryGeneratorInfo.h:100
Generated on Wed Mar 20 2024 12:10:22 for REST-for-Physics by doxygen 1.9.4