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 
16 using namespace std;
17 using namespace TRestGeant4PrimaryGeneratorTypes;
18 
19 string 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 
40 SpatialGeneratorTypes 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 
67 string 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 
88 SpatialGeneratorShapes 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 
116 string 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 
138 EnergyDistributionTypes 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 
166 string 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 
180 EnergyDistributionFormulas 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 
197 TF1 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 
224 string 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 
246 AngularDistributionTypes 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 
274 string 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 
288 AngularDistributionFormulas 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 
304 TF1 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 
338 string 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 
350 EnergyAndAngularDistributionFormulas
351 TRestGeant4PrimaryGeneratorTypes::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 
364 TF2 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 
402 void 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 }
Generated on Thu Nov 7 2024 11:22:21 for REST-for-Physics by doxygen 1.9.1