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