Skip to content

Commit af24dd5

Browse files
lobislobis
authored
Merge pull request #124 from rest-for-physics/lobis-cosmic-revamp
cosmic sampler
2 parents 2bac47b + db81050 commit af24dd5

7 files changed

Lines changed: 59 additions & 17 deletions

‎inc/TRestGeant4ParticleSourceCosmics.h‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,8 @@ class TRestGeant4ParticleSourceCosmics : public TRestGeant4ParticleSource {
1313
std::map<std::string, double> fParticleWeights;
1414

1515
std::map<std::string, TH2D*> fHistograms;
16+
std::map<std::string, TH2D*> fHistogramsTransformed;
17+
1618
static std::mutex fMutex;
1719
static std::unique_ptr<TRandom3> fRandom;
1820

@@ -27,7 +29,9 @@ class TRestGeant4ParticleSourceCosmics : public TRestGeant4ParticleSource {
2729

2830
const char* GetName() const override { return "TRestGeant4ParticleSourceCosmics"; }
2931

30-
ClassDefOverride(TRestGeant4ParticleSourceCosmics, 1);
32+
std::map<std::string, TH2D*> GetHistogramsTransformed() const { return fHistogramsTransformed; }
33+
34+
ClassDefOverride(TRestGeant4ParticleSourceCosmics, 2);
3135
};
3236

3337
#endif // REST_TRESTGEANT4PARTICLESOURCECOSMICS_H

‎inc/TRestGeant4PrimaryGeneratorInfo.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ SpatialGeneratorShapes StringToSpatialGeneratorShapes(const std::string&);
3737

3838
enum class EnergyDistributionTypes {
3939
TH1D,
40+
TH2D,
4041
FORMULA,
4142
FORMULA2,
4243
MONO,
@@ -58,6 +59,7 @@ TF1 EnergyDistributionFormulasToRootFormula(const EnergyDistributionFormulas&);
5859

5960
enum class AngularDistributionTypes {
6061
TH1D,
62+
TH2D,
6163
FORMULA,
6264
FORMULA2,
6365
ISOTROPIC,

‎src/TRestGeant4AnalysisProcess.cxx‎

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -514,10 +514,15 @@ TRestEvent* TRestGeant4AnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
514514
}
515515
SetObservableValue("subEventPrimaryParticleName", subEventPrimaryParticleName);
516516

517-
TVector3 v(xDirection, yDirection, zDirection);
518-
SetObservableValue("thetaPrimary", v.Theta());
519-
SetObservableValue("phiPrimary", v.Phi());
517+
TVector3 primaryDirection(xDirection, yDirection, zDirection);
518+
primaryDirection = primaryDirection.Unit();
519+
SetObservableValue("thetaPrimary", primaryDirection.Theta());
520+
SetObservableValue("phiPrimary", primaryDirection.Phi());
520521
SetObservableValue("zenithYDegrees", TMath::ACos(-1 * yDirection) * TMath::RadToDeg());
522+
TVector3 sourceDirection = fG4Metadata->GetParticleSource(0)->GetDirection();
523+
sourceDirection = sourceDirection.Unit();
524+
SetObservableValue("zenithSourceDegrees",
525+
TMath::ACos(primaryDirection.Dot(sourceDirection)) * TMath::RadToDeg());
521526

522527
Double_t energyPrimary = fOutputG4Event->GetPrimaryEventEnergy(0);
523528
SetObservableValue("energyPrimary", energyPrimary);

‎src/TRestGeant4Metadata.cxx‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1110,8 +1110,10 @@ void TRestGeant4Metadata::ReadParticleSource(TRestGeant4ParticleSource* source,
11101110
}
11111111
source->SetAngularDistributionType(GetParameter(
11121112
"type", angularDefinition, AngularDistributionTypesToString(AngularDistributionTypes::FLUX)));
1113-
if (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1114-
AngularDistributionTypes::TH1D) {
1113+
if ((StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1114+
AngularDistributionTypes::TH1D) ||
1115+
(StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1116+
AngularDistributionTypes::TH2D)) {
11151117
source->SetAngularDistributionFilename(SearchFile(GetParameter("file", angularDefinition)));
11161118
auto name = GetParameter("name", angularDefinition);
11171119
if (name == "NO_SUCH_PARA") {
@@ -1164,8 +1166,10 @@ void TRestGeant4Metadata::ReadParticleSource(TRestGeant4ParticleSource* source,
11641166
}
11651167
source->SetEnergyDistributionType(GetParameter(
11661168
"type", energyDefinition, EnergyDistributionTypesToString(EnergyDistributionTypes::MONO)));
1167-
if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1168-
EnergyDistributionTypes::TH1D) {
1169+
if ((StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1170+
EnergyDistributionTypes::TH1D) ||
1171+
(StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1172+
EnergyDistributionTypes::TH2D)) {
11691173
source->SetEnergyDistributionFilename(SearchFile(GetParameter("file", energyDefinition)));
11701174
auto name = GetParameter("name", energyDefinition);
11711175
if (name == "NO_SUCH_PARA") {

‎src/TRestGeant4ParticleSource.cxx‎

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -47,8 +47,10 @@ void TRestGeant4ParticleSource::PrintMetadata() {
4747
RESTMetadata << "Particle charge: " << GetParticleCharge() << RESTendl;
4848
}
4949
RESTMetadata << "Angular distribution type: " << GetAngularDistributionType() << RESTendl;
50-
if (StringToAngularDistributionTypes(GetAngularDistributionType().Data()) ==
51-
AngularDistributionTypes::TH1D) {
50+
if ((StringToAngularDistributionTypes(GetAngularDistributionType().Data()) ==
51+
AngularDistributionTypes::TH1D) ||
52+
(StringToAngularDistributionTypes(GetAngularDistributionType().Data()) ==
53+
AngularDistributionTypes::TH2D)) {
5254
RESTMetadata << "Angular distribution filename: "
5355
<< TRestTools::GetPureFileName((string)GetAngularDistributionFilename()) << RESTendl;
5456
RESTMetadata << "Angular distribution histogram name: " << GetAngularDistributionNameInFile()
@@ -68,18 +70,20 @@ void TRestGeant4ParticleSource::PrintMetadata() {
6870
}
6971
}
7072
if ((StringToAngularDistributionTypes(GetAngularDistributionType().Data()) ==
71-
AngularDistributionTypes::FORMULA) ||
73+
AngularDistributionTypes::TH1D) ||
7274
StringToAngularDistributionTypes(GetAngularDistributionType().Data()) ==
73-
AngularDistributionTypes::FORMULA2) {
75+
AngularDistributionTypes::TH2D) {
7476
RESTMetadata << "Angular distribution range (deg): ("
7577
<< GetAngularDistributionRangeMin() * TMath::RadToDeg() << ", "
7678
<< GetAngularDistributionRangeMax() * TMath::RadToDeg() << ")" << RESTendl;
7779
RESTMetadata << "Angular distribution random sampling grid size: "
7880
<< GetAngularDistributionFormulaNPoints() << RESTendl;
7981
}
8082
RESTMetadata << "Energy distribution type: " << GetEnergyDistributionType() << RESTendl;
81-
if (StringToEnergyDistributionTypes(GetEnergyDistributionType().Data()) ==
82-
EnergyDistributionTypes::TH1D) {
83+
if ((StringToEnergyDistributionTypes(GetEnergyDistributionType().Data()) ==
84+
EnergyDistributionTypes::TH1D) ||
85+
StringToEnergyDistributionTypes(GetEnergyDistributionType().Data()) ==
86+
EnergyDistributionTypes::TH2D) {
8387
RESTMetadata << "Energy distribution filename: "
8488
<< TRestTools::GetPureFileName((string)GetEnergyDistributionFilename()) << RESTendl;
8589
RESTMetadata << "Energy distribution histogram name: " << GetEnergyDistributionNameInFile()

‎src/TRestGeant4ParticleSourceCosmics.cxx‎

Lines changed: 16 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,8 @@ void TRestGeant4ParticleSourceCosmics::InitFromConfigFile() {
2121
fFilename = GetParameter("filename");
2222
const auto particles = GetParameter("particles", "neutron,proton,gamma,electron,muon");
2323

24+
fDirection = Get3DVectorParameterWithUnits("direction", {0, -1, 0});
25+
2426
std::istringstream iss(particles);
2527
std::string name;
2628
while (std::getline(iss, name, ',')) {
@@ -84,8 +86,20 @@ void TRestGeant4ParticleSourceCosmics::Update() {
8486
}
8587
}
8688

87-
auto& hist = fHistograms.at(particleName);
88-
89+
if (fHistogramsTransformed.find(particleName) == fHistogramsTransformed.end()) {
90+
auto histOriginal = fHistograms.at(particleName);
91+
TH2D* hist = new TH2D(*histOriginal);
92+
// same as original but Y axis is multiplied by 1/cos(zenith). Integral should be the same.
93+
for (int i = 1; i <= hist->GetNbinsX(); i++) {
94+
for (int j = 1; j <= hist->GetNbinsY(); j++) {
95+
const double zenith = hist->GetYaxis()->GetBinCenter(j);
96+
const double value = hist->GetBinContent(i, j) / TMath::Cos(zenith * TMath::DegToRad());
97+
hist->SetBinContent(i, j, value);
98+
}
99+
}
100+
fHistogramsTransformed[particleName] = hist;
101+
}
102+
auto hist = fHistogramsTransformed.at(particleName);
89103
double energy, zenith;
90104
hist->GetRandom2(energy, zenith, fRandom.get());
91105

‎src/TRestGeant4PrimaryGeneratorInfo.cxx‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,6 @@
99
#include <TF2.h>
1010
#include <TMath.h>
1111
#include <TRestStringOutput.h>
12-
#include <tinyxml.h>
1312

1413
#include <iostream>
1514

@@ -118,6 +117,8 @@ string TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypesToString(
118117
switch (type) {
119118
case EnergyDistributionTypes::TH1D:
120119
return "TH1D";
120+
case EnergyDistributionTypes::TH2D:
121+
return "TH2D";
121122
case EnergyDistributionTypes::FORMULA:
122123
return "Formula";
123124
case EnergyDistributionTypes::FORMULA2:
@@ -140,6 +141,9 @@ EnergyDistributionTypes TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistribu
140141
if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::TH1D),
141142
TString::ECaseCompare::kIgnoreCase)) {
142143
return EnergyDistributionTypes::TH1D;
144+
} else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::TH2D),
145+
TString::ECaseCompare::kIgnoreCase)) {
146+
return EnergyDistributionTypes::TH2D;
143147
} else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::FORMULA),
144148
TString::ECaseCompare::kIgnoreCase)) {
145149
return EnergyDistributionTypes::FORMULA;
@@ -226,6 +230,8 @@ string TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypesToString(
226230
switch (type) {
227231
case AngularDistributionTypes::TH1D:
228232
return "TH1D";
233+
case AngularDistributionTypes::TH2D:
234+
return "TH2D";
229235
case AngularDistributionTypes::FORMULA:
230236
return "Formula";
231237
case AngularDistributionTypes::FORMULA2:
@@ -248,6 +254,9 @@ AngularDistributionTypes TRestGeant4PrimaryGeneratorTypes::StringToAngularDistri
248254
if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::TH1D),
249255
TString::ECaseCompare::kIgnoreCase)) {
250256
return AngularDistributionTypes::TH1D;
257+
} else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::TH2D),
258+
TString::ECaseCompare::kIgnoreCase)) {
259+
return AngularDistributionTypes::TH2D;
251260
} else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::FORMULA),
252261
TString::ECaseCompare::kIgnoreCase)) {
253262
return AngularDistributionTypes::FORMULA;

0 commit comments

Comments
 (0)