Skip to content

Commit

Permalink
Merge branch 'master' into lobis-cosmic-revamp
Browse files Browse the repository at this point in the history
  • Loading branch information
lobis authored Apr 9, 2024
2 parents 1330982 + f1c00d2 commit 2ae86b3
Show file tree
Hide file tree
Showing 8 changed files with 186 additions and 39 deletions.
4 changes: 2 additions & 2 deletions .github/workflows/validation.yml
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,7 @@ jobs:
source ${{ env.REST_INSTALL_PATH }}/thisREST.sh
cd ${{ env.REST_INSTALL_PATH }}/examples/restG4/07.FullChainDecay/
restG4 fullChain.rml -o Run00001_U238_FullChainDecay.root
restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 15)'
restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 16)'
restG4 singleDecay.rml -o Run00002_U238_SingleChainDecay.root
restRoot -b -q Validate.C'("Run00002_U238_SingleChainDecay.root", 1)'
export REST_ISOTOPE=Be7
Expand Down Expand Up @@ -538,7 +538,7 @@ jobs:
source ${{ env.REST_INSTALL_PATH }}/thisREST.sh
cd ${{ env.REST_INSTALL_PATH }}/examples/restG4/07.FullChainDecay/
restG4 fullChain.rml -o Run00001_U238_FullChainDecay.root
restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 15)'
restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 16)'
restG4 singleDecay.rml -o Run00002_U238_SingleChainDecay.root
restRoot -b -q Validate.C'("Run00002_U238_SingleChainDecay.root", 1)'
export REST_ISOTOPE=Be7
Expand Down
2 changes: 1 addition & 1 deletion .gitlab-ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@ Load REST Libraries:
- cd ${CI_PROJECT_DIR}/install/examples/restG4/07.FullChainDecay/
- restG4 fullChain.rml
- restG4 singleDecay.rml
- restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 15)'
- restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 16)'
- restRoot -b -q Validate.C'("Run00002_U238_SingleChainDecay.root", 1)'
- export REST_ISOTOPE="Be7"
- restG4 singleDecay.rml
Expand Down
2 changes: 1 addition & 1 deletion examples/06.IonRecoils/Validate.C
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ Int_t Validate(const char* filename) {
const bool isReferenceGeant4Version = geant4Metadata->GetGeant4Version() == "10.4.3";

double averageNumberOfHits = 0;
const double averageNumberOfHitsRef = (!isReferenceGeant4Version) ? 14417.5 : 11223.0;
const double averageNumberOfHitsRef = (!isReferenceGeant4Version) ? 11131.5 : 10369.5;
const double tolerance = 0.001;

for (int i = 0; i < run.GetEntries(); i++) {
Expand Down
10 changes: 5 additions & 5 deletions examples/08.Alphas/Validate.C
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ Int_t Validate() {
run.GetEntry(100);

Double_t thetaSample = analysisTree->GetObservableValue<Double_t>("g4Ana_thetaPrimary");
constexpr Double_t thetaSampleRef = 2.85884;
constexpr Double_t thetaSampleRef = 2.29387;

Double_t phiSample = analysisTree->GetObservableValue<Double_t>("g4Ana_phiPrimary");
constexpr Double_t phiSampleRef = -1.59582;
Expand All @@ -24,24 +24,24 @@ Int_t Validate() {
constexpr Double_t thetaAverageRef = 1.5767;

Double_t phiAverage = analysisTree->GetObservableAverage("g4Ana_phiPrimary");
constexpr Double_t phiAverageRef = 0.0719113;
constexpr Double_t phiAverageRef = -0.0169888;

cout << "entry 100, theta: " << thetaSample << ", phi: " << phiSample << endl;
cout << "average theta: " << thetaAverage << ", average phi: " << phiAverage << endl;

if (h1->Integral() != 1984) {
if (h1->Integral() != 1926) {
cout << "Wrong number of alphas produced at (h1) 1MeV_5um!" << endl;
cout << "Histogram contains : " << h1->Integral() << endl;
return 10;
}

if (h2->Integral() != 7432) {
if (h2->Integral() != 7487) {
cout << "Wrong number of alphas produced at (h2) 5MeV_5um!" << endl;
cout << "Histogram contains : " << h2->Integral() << endl;
return 20;
}

if (h3->Integral() != 9430) {
if (h3->Integral() != 9432) {
cout << "Wrong number of alphas produced at (h3) 5MeV_1um!" << endl;
cout << "Histogram contains : " << h3->Integral() << endl;
return 30;
Expand Down
43 changes: 43 additions & 0 deletions examples/12.Generators/IsotropicWithRange.rml
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
<?xml version="1.0" encoding="UTF-8" standalone="no" ?>

<restG4>
<globals>
<parameter name="mainDataPath" value=""/>
<parameter name="verboseLevel" value="essential"/>
</globals>

<TRestRun>
<parameter name="experimentName" value="IsotropicWithRange"/>
<parameter name="runType" value="simulation"/>
<parameter name="runNumber" value="1"/>
<parameter name="runTag" value="IsotropicWithRange"/>
<parameter name="outputFileName" value="IsotropicWithRange.root"/>
<parameter name="runDescription" value=""/>
<parameter name="user" value="${USER}"/>
<parameter name="overwrite" value="off"/>
<parameter name="readOnly" value="false"/>
</TRestRun>

<TRestGeant4Metadata>
<parameter name="gdmlFile" value="geometry/setup.gdml"/>
<parameter name="seed" value="171981"/>
<parameter name="nEvents" value="1000000"/>
<parameter name="saveAllEvents" value="true"/>
<generator type="surface" shape="circle"
position="(0,100,0)mm" size="(400,0,0)mm"
rotationAngle="90deg" rotationAxis="(1,0,0)">
<source particle="geantino">
<angular type="isotropic" direction="(0,-1,0)" coneHalfAngle="15deg"/>
<energy type="mono" energy="10.0" units="MeV"/>
</source>
</generator>

<detector>
<parameter name="energyRange" value="(0,10)" units="GeV"/>
<volume name="gasVolume" sensitive="true" maxStepSize="1mm"/>
</detector>

</TRestGeant4Metadata>

<TRestGeant4PhysicsLists name="default" title="First physics list implementation"/>
</restG4>
94 changes: 94 additions & 0 deletions examples/12.Generators/ValidateIsotropicWithRange.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
#include <TRestGeant4Event.h>

Int_t ValidateIsotropicWithRange(const char* filename) {
cout << "Starting validation for '" << filename << "'" << endl;

TRestRun run(filename);

cout << "Run entries: " << run.GetEntries() << endl;
if (run.GetEntries() != 1000000) {
cout << "Bad number of entries: " << run.GetEntries() << endl;
return 2;
}

auto metadata = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata");
TRestGeant4Event* event = run.GetInputEvent<TRestGeant4Event>();
TVector3 sourceDirection = metadata->GetParticleSource()->GetDirection();
cout << "original source direction: (" << sourceDirection.X() << ", " << sourceDirection.Y() << ", "
<< sourceDirection.Z() << ")" << endl;

constexpr double tolerance = 0.1;

double thetaAverage = 0, thetaMin = TMath::Infinity(), thetaMax = 0;
constexpr double thetaAverageRef = 10.0, thetaMinRef = 0.0, thetaMaxRef = 15.0;

double energyPrimaryAverage = 0, energyPrimaryMin = TMath::Infinity(), energyPrimaryMax = 0;
constexpr double energyPrimaryAverageRef = 10000.0, energyPrimaryMinRef = 10000.0,
energyPrimaryMaxRef = 10000.0;

for (int i = 0; i < run.GetEntries(); i++) {
run.GetEntry(i);
TVector3 direction = event->GetPrimaryEventDirection();
const auto theta = sourceDirection.Angle(direction) * TMath::RadToDeg();

thetaAverage += theta / run.GetEntries();
if (theta < thetaMin) {
thetaMin = theta;
}
if (theta > thetaMax) {
thetaMax = theta;
}

Double_t energy = event->GetPrimaryEventEnergy();
energyPrimaryAverage += energy / run.GetEntries();
if (energy < energyPrimaryMin) {
energyPrimaryMin = energy;
}
if (energy > energyPrimaryMax) {
energyPrimaryMax = energy;
}
}

cout << "Average theta (deg): " << thetaAverage << endl;
cout << "Minimum theta (deg): " << thetaMin << endl;
cout << "Maximum theta (deg): " << thetaMax << endl;

cout << "Average energy (keV): " << energyPrimaryAverage << endl;
cout << "Minimum energy (keV): " << energyPrimaryMin << endl;
cout << "Maximum energy (keV): " << energyPrimaryMax << endl;

if (TMath::Abs(thetaAverage - thetaAverageRef) / thetaAverageRef > tolerance) {
cout << "The average theta of the distribution is wrong!" << endl;
cout << "thetaAverage (deg): " << thetaAverage << endl;
return 3;
}
if (TMath::Abs(thetaMin - thetaMinRef) / (thetaMinRef + 1) > tolerance) {
cout << "The minimum theta of the distribution is wrong!" << endl;
cout << "thetaMin (deg): " << thetaMin << endl;
return 4;
}
if (TMath::Abs(thetaMax - thetaMaxRef) / thetaMaxRef > tolerance) {
cout << "The maximum theta of the distribution is wrong!" << endl;
cout << "thetaMax (deg): " << thetaMax << endl;
return 5;
}

if (TMath::Abs(energyPrimaryAverage - energyPrimaryAverageRef) / energyPrimaryAverageRef > tolerance) {
cout << "The average energy of the distribution is wrong!" << endl;
cout << "energyPrimaryAverage (keV): " << energyPrimaryAverage << endl;
return 6;
}
if (TMath::Abs(energyPrimaryMin - energyPrimaryMinRef) / energyPrimaryMinRef > tolerance) {
cout << "The minimum energy of the distribution is wrong!" << endl;
cout << "energyPrimaryMin (keV): " << energyPrimaryMin << endl;
return 7;
}
if (TMath::Abs(energyPrimaryMax - energyPrimaryMaxRef) / energyPrimaryMaxRef > tolerance) {
cout << "The maximum energy of the distribution is wrong!" << endl;
cout << "energyPrimaryMax (keV): " << energyPrimaryMax << endl;
return 8;
}

cout << "All tests passed! [\033[32mOK\033[0m]\n";
return 0;
}
44 changes: 15 additions & 29 deletions src/PrimaryGeneratorAction.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -421,8 +421,8 @@ void PrimaryGeneratorAction::SetParticleDirection(Int_t particleSourceIndex,
TRestGeant4Metadata* restG4Metadata = simulationManager->GetRestMetadata();
TRestGeant4ParticleSource* source = restG4Metadata->GetParticleSource(0);

G4ThreeVector direction = {source->GetDirection().X(), source->GetDirection().Y(),
source->GetDirection().Z()};
const auto sourceDirection = source->GetDirection();
G4ThreeVector direction = {sourceDirection.X(), sourceDirection.Y(), sourceDirection.Z()};

const string angularDistTypeName = source->GetAngularDistributionType().Data();
const auto angularDistTypeEnum = StringToAngularDistributionTypes(angularDistTypeName);
Expand All @@ -431,20 +431,16 @@ void PrimaryGeneratorAction::SetParticleDirection(Int_t particleSourceIndex,
cout << "DEBUG: Angular distribution: " << angularDistTypeName << endl;
}

// generator type
/* Apparently not used
const auto& primaryGeneratorInfo = restG4Metadata->GetGeant4PrimaryGeneratorInfo();
const string& spatialGeneratorTypeName = primaryGeneratorInfo.GetSpatialGeneratorType().Data();
const auto spatialGeneratorTypeEnum = StringToSpatialGeneratorTypes(spatialGeneratorTypeName);
const string& spatialGeneratorShapeName = primaryGeneratorInfo.GetSpatialGeneratorShape().Data();
const auto spatialGeneratorShapeEnum = StringToSpatialGeneratorShapes(spatialGeneratorShapeName);
*/

if (angularDistTypeEnum == AngularDistributionTypes::ISOTROPIC) {
direction = GetIsotropicVector();
if (source->GetAngularDistributionIsotropicConeHalfAngle() > 0) {
const auto originalDirection = direction;
do {
direction = GetIsotropicVector();
} while (originalDirection.angle(direction) >
source->GetAngularDistributionIsotropicConeHalfAngle());
} else {
direction = GetIsotropicVector();
}
} else if (angularDistTypeEnum == AngularDistributionTypes::TH1D) {
Double_t angle = 0;
Double_t value = G4UniformRand() * fAngularDistributionHistogram->Integral();
Expand Down Expand Up @@ -656,20 +652,10 @@ void PrimaryGeneratorAction::SetParticlePosition() {
}

G4ThreeVector PrimaryGeneratorAction::GetIsotropicVector() const {
G4double a, b, c;
G4double n;
do {
a = (G4UniformRand() - 0.5) / 0.5;
b = (G4UniformRand() - 0.5) / 0.5;
c = (G4UniformRand() - 0.5) / 0.5;
n = a * a + b * b + c * c;
} while (n > 1 || n == 0.0);

n = sqrt(n);
a /= n;
b /= n;
c /= n;
return {a, b, c};
double phi = 2 * M_PI * G4UniformRand();
double theta = TMath::ACos(1 - 2 * G4UniformRand());

return {TMath::Sin(theta) * TMath::Cos(phi), TMath::Sin(theta) * TMath::Sin(phi), TMath::Cos(theta)};
}

void PrimaryGeneratorAction::GenPositionOnGDMLVolume(double& x, double& y, double& z) {
Expand Down
26 changes: 25 additions & 1 deletion test/src/examples.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ TEST(restG4, Example_07_Decay_FullChain) {
gROOT->ProcessLine(TString::Format(".L %s", macro.Data())); // Load macro
int error = 0;
const int result =
gROOT->ProcessLine(TString::Format("Validate(\"%s\", %d)", options.outputFile.c_str(), 15), &error);
gROOT->ProcessLine(TString::Format("Validate(\"%s\", %d)", options.outputFile.c_str(), 16), &error);
EXPECT_EQ(error, 0);
EXPECT_EQ(result, 0);

Expand Down Expand Up @@ -417,6 +417,30 @@ TEST(restG4, Example_12_Generators_EnergyAndAngularCorrelated) {
fs::current_path(originalPath);
}

TEST(restG4, Example_12_Generators_IsotropicWithRange) {
const auto originalPath = fs::current_path();
const auto thisExamplePath = examplesPath / "12.Generators";
fs::current_path(thisExamplePath);

CommandLineOptions::Options options;
options.rmlFile = "IsotropicWithRange.rml";
options.outputFile = thisExamplePath / "IsotropicWithRange.root";

Application app;
app.Run(options);

// Run validation macro
const TString macro(thisExamplePath / "ValidateIsotropicWithRange.C");
gROOT->ProcessLine(TString::Format(".L %s", macro.Data())); // Load macro
int error = 0;
const int result = gROOT->ProcessLine(
TString::Format("ValidateIsotropicWithRange(\"%s\")", options.outputFile.c_str()), &error);
EXPECT_EQ(error, 0);
EXPECT_EQ(result, 0);

fs::current_path(originalPath);
}

TEST(restG4, Example_13_IAXO_Neutrons) {
// cd into example
const auto originalPath = fs::current_path();
Expand Down

0 comments on commit 2ae86b3

Please sign in to comment.