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