update cosmic generator example

examples/12.Generators/CosmicGenerator.rml

@@ -7,11 +7,11 @@
     </globals>
 
     <TRestRun>
-        <parameter name="experimentName" value="CosmicMuons"/>
+        <parameter name="experimentName" value="CosmicGenerator"/>
         <parameter name="runType" value="simulation"/>
         <parameter name="runNumber" value="1"/>
-        <parameter name="runTag" value="CosmicMuons"/>
-        <parameter name="outputFileName" value="CosmicMuons.root"/>
+        <parameter name="runTag" value="CosmicGenerator"/>
+        <parameter name="outputFileName" value="CosmicGenerator.root"/>
         <parameter name="runDescription" value=""/>
         <parameter name="user" value="${USER}"/>
         <parameter name="overwrite" value="off"/>
@@ -21,17 +21,16 @@
     <TRestGeant4Metadata>
         <parameter name="gdmlFile" value="geometry/setup.gdml"/>
         <parameter name="seed" value="137"/>
-        <parameter name="nRequestedEntries" value="100000"/>
-        <parameter name="saveAllEvents" value="false"/>
+        <parameter name="nRequestedEntries" value="1000000"/>
+        <parameter name="saveAllEvents" value="true"/>
         <generator type="cosmic">
             <source particle="geantino">
-                <energy type="formula2" name="CosmicMuons" nPoints="1200"/>
-                <angular type="formula2" direction="(0,-1,0)" nPoints="300"/>
+                <energy type="formula2" name="CosmicMuons" nPoints="2000"/>
+                <angular type="formula2" direction="(0.1,-0.5,1.25)" nPoints="500"/>
             </source>
         </generator>
 
         <detector>
-            <parameter name="energyRange" value="(0,10)" units="GeV"/>
             <volume name="gasVolume" sensitive="true" maxStepSize="1mm"/>
         </detector>
 

examples/12.Generators/CosmicMuonsEnergyAngularCorrelated.rml

@@ -7,11 +7,11 @@
     </globals>
 
     <TRestRun>
-        <parameter name="experimentName" value="CosmicMuons"/>
+        <parameter name="experimentName" value="CosmicMuonsEnergyAngularCorrelated"/>
         <parameter name="runType" value="simulation"/>
         <parameter name="runNumber" value="1"/>
-        <parameter name="runTag" value="CosmicMuons"/>
-        <parameter name="outputFileName" value="CosmicMuons.root"/>
+        <parameter name="runTag" value="CosmicMuonsEnergyAngularCorrelated"/>
+        <parameter name="outputFileName" value="CosmicMuonsEnergyAngularCorrelated.root"/>
         <parameter name="runDescription" value=""/>
         <parameter name="user" value="${USER}"/>
         <parameter name="overwrite" value="off"/>
@@ -21,14 +21,14 @@
     <TRestGeant4Metadata>
         <parameter name="gdmlFile" value="geometry/setup.gdml"/>
         <parameter name="seed" value="171981"/>
-        <parameter name="nEvents" value="100000"/>
+        <parameter name="nEvents" value="200000"/>
         <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">
                 <!-- Use a TF2 to sample energy and angle at once (correlation) -->
-                <energy type="formula2" name="CosmicMuons" range="(100,150)MeV" nPoints="1200"/>
+                <energy type="formula2" name="CosmicMuons" range="(100,850)MeV" nPoints="1200"/>
                 <angular type="formula2" direction="(0,-1,0)"
                          range="(50,20)deg" nPoints="300"/> <!-- name is not required as it is redundant -->
             </source>

examples/12.Generators/ValidateCosmicGenerator.C

@@ -6,31 +6,42 @@ Int_t ValidateCosmicGenerator(const char* filename) {
     TRestRun run(filename);
 
     cout << "Run entries: " << run.GetEntries() << endl;
-    if (run.GetEntries() != 100000) {
+    if (run.GetEntries() != 1E6) {
         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;
 
+    TVector3 sourceDirection = metadata->GetParticleSource()->GetDirection();
+    if (TVector3(0.1, -0.5, 1.25).Angle(sourceDirection) > 0.01) {
+        cout << "Bad source direction: (" << sourceDirection.X() << ", " << sourceDirection.Y() << ", "
+             << sourceDirection.Z() << endl;
+        return 3;
+    }
+
     double thetaAverage = 0, thetaMin = TMath::Infinity(), thetaMax = 0;
-    constexpr double thetaAverageRef = 45.39, thetaMinRef = 0.32, thetaMaxRef = 89.99;
+    constexpr double thetaAverageRef = 37.8595, thetaMinRef = 0.0110538, thetaMaxRef = 89.9982;
 
     double energyPrimaryAverage = 0, energyPrimaryMin = TMath::Infinity(), energyPrimaryMax = 0;
-    constexpr double energyPrimaryAverageRef = 8590750, energyPrimaryMinRef = 145.8,
-                     energyPrimaryMaxRef = 2810890000;
+    constexpr double energyPrimaryAverageRef = 7.89997e+06, energyPrimaryMinRef = 200008,
+                     energyPrimaryMaxRef = 4.63073e+09;
+
+    double primaryRadiusRef = 173.20508;
 
     TH1D thetaHist("thetaHist", "Theta angle from source direction", 100, 0, TMath::Pi() * TMath::RadToDeg());
     TH1D energyHist("energyHist", "Primary muon energy", 200, 0, 1E9);
     for (int i = 0; i < run.GetEntries(); i++) {
         run.GetEntry(i);
         TVector3 direction = event->GetPrimaryEventDirection();
+        TVector3 position = event->GetPrimaryEventOrigin();
+        if (TMath::Abs(position.Mag() - primaryRadiusRef) / primaryRadiusRef > tolerance) {
+            cout << "Bad primary position: (" << position.X() << ", " << position.Y() << ", " << position.Z()
+                 << ") with radius " << position.Mag() << endl;
+            return 4;
+        }
         const auto theta = sourceDirection.Angle(direction) * TMath::RadToDeg();
         thetaHist.Fill(theta);
 
@@ -67,52 +78,60 @@ Int_t ValidateCosmicGenerator(const char* filename) {
     if (TMath::Abs(thetaAverage - thetaAverageRef) / thetaAverageRef > tolerance) {
         cout << "The average theta of the distribution is wrong!" << endl;
         cout << "thetaAverage (deg): " << thetaAverage << endl;
-        return 3;
+        return 5;
     }
     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;
+        return 6;
     }
     if (TMath::Abs(thetaMax - thetaMaxRef) / thetaMaxRef > tolerance) {
         cout << "The maximum theta of the distribution is wrong!" << endl;
         cout << "thetaMax (deg): " << thetaMax << endl;
-        return 5;
+        return 7;
     }
 
     if (TMath::Abs(energyPrimaryAverage - energyPrimaryAverageRef) / energyPrimaryAverageRef > tolerance) {
         cout << "The average energy of the distribution is wrong!" << endl;
         cout << "energyPrimaryAverage (keV): " << energyPrimaryAverage << endl;
-        return 6;
+        return 8;
     }
     if (TMath::Abs(energyPrimaryMin - energyPrimaryMinRef) / energyPrimaryMinRef > tolerance) {
         cout << "The minimum energy of the distribution is wrong!" << endl;
         cout << "energyPrimaryMin (keV): " << energyPrimaryMin << endl;
-        return 7;
+        return 9;
     }
     if (TMath::Abs(energyPrimaryMax - energyPrimaryMaxRef) / energyPrimaryMaxRef > tolerance) {
         cout << "The maximum energy of the distribution is wrong!" << endl;
         cout << "energyPrimaryMax (keV): " << energyPrimaryMax << endl;
-        return 8;
+        return 10;
     }
 
-    if (metadata->GetNumberOfEvents() != 885919) {
+    if (metadata->GetNumberOfEvents() != 1E6) {
         cout << "wrong number of events: " << metadata->GetNumberOfEvents() << endl;
-        return 9;
+        return 11;
     }
 
     const auto surfaceTermRef = 942.4778;
-    if ((metadata->GetGeant4PrimaryGeneratorInfo().GetSpatialGeneratorCosmicSurfaceTerm() - surfaceTermRef) /
+    if (abs(metadata->GetGeant4PrimaryGeneratorInfo().GetSpatialGeneratorCosmicSurfaceTermCm2() -
+            surfaceTermRef) /
             surfaceTermRef >
         tolerance) {
         cout << "wrong cosmic surface term: "
-             << metadata->GetGeant4PrimaryGeneratorInfo().GetSpatialGeneratorCosmicSurfaceTerm() << endl;
-        return 10;
+             << metadata->GetGeant4PrimaryGeneratorInfo().GetSpatialGeneratorCosmicSurfaceTermCm2() << endl;
+        return 12;
+    }
+
+    const auto cosmicFluxRef = 0.108215;
+    if (abs(metadata->GetCosmicFluxInCountsPerCm2PerSecond() - cosmicFluxRef) / cosmicFluxRef > tolerance) {
+        cout << "wrong cosmic flux: " << metadata->GetCosmicFluxInCountsPerCm2PerSecond() << endl;
+        return 13;
     }
-    const auto simulationTimeRef = 71193.790;
-    if ((metadata->GetEquivalentSimulatedTime() - simulationTimeRef) / simulationTimeRef > tolerance) {
+
+    const auto simulationTimeRef = 9804.87;
+    if (abs(metadata->GetEquivalentSimulatedTime() - simulationTimeRef) / simulationTimeRef > tolerance) {
         cout << "wrong equivalent simulation time: " << metadata->GetEquivalentSimulatedTime() << endl;
-        return 11;
+        return 14;
     }
 
     cout << "All tests passed! [\033[32mOK\033[0m]\n";

examples/12.Generators/ValidateCosmicMuonsEnergyAngularCorrelated.C

@@ -6,7 +6,7 @@ Int_t ValidateCosmicMuonsEnergyAngularCorrelated(const char* filename) {
     TRestRun run(filename);
 
     cout << "Run entries: " << run.GetEntries() << endl;
-    if (run.GetEntries() != 100000) {
+    if (run.GetEntries() != 200000) {
         cout << "Bad number of entries: " << run.GetEntries() << endl;
         return 2;
     }
@@ -20,11 +20,11 @@ Int_t ValidateCosmicMuonsEnergyAngularCorrelated(const char* filename) {
     constexpr double tolerance = 0.1;
 
     double thetaAverage = 0, thetaMin = TMath::Infinity(), thetaMax = 0;
-    constexpr double thetaAverageRef = 33.8542, thetaMinRef = 20.0, thetaMaxRef = 50.0;
+    constexpr double thetaAverageRef = 34.0558, thetaMinRef = 20., thetaMaxRef = 50.;
 
     double energyPrimaryAverage = 0, energyPrimaryMin = TMath::Infinity(), energyPrimaryMax = 0;
-    constexpr double energyPrimaryAverageRef = 124804.0, energyPrimaryMinRef = 100000,
-                     energyPrimaryMaxRef = 150000.0;
+    constexpr double energyPrimaryAverageRef = 507020., energyPrimaryMinRef = 200000,
+                     energyPrimaryMaxRef = 850000;
 
     TH1D thetaHist("thetaHist", "Theta angle from source direction", 100, 0, TMath::Pi() * TMath::RadToDeg());
     TH1D energyHist("energyHist", "Primary muon energy", 200, 0, 1E9);

src/PrimaryGeneratorAction.cxx

@@ -23,22 +23,18 @@ using namespace std;
 using namespace TRestGeant4PrimaryGeneratorTypes;
 
 G4ThreeVector ComputeCosmicPosition(const G4ThreeVector& direction, double radius) {
-    // Angles in radians
-    const auto directionRoot = TVector3(direction.x(), direction.y(), direction.z());
-    const auto theta = directionRoot.Angle(TVector3(0, -1, 0));
-    const auto phi = directionRoot.Phi();
-    // Get random point in a disk
+    // Get random point in a disk with 'direction' as normal
     const double u1 = G4UniformRand(), u2 = G4UniformRand();
     const auto positionInDisk =
-        G4ThreeVector(sqrt(u1) * cos(2. * M_PI * u2), 0, sqrt(u1) * sin(2. * M_PI * u2))
-            .rotateX(theta)
-            .rotateY(phi) *
+        G4ThreeVector(sqrt(u1) * cos(2. * M_PI * u2), sqrt(u1) * sin(2. * M_PI * u2), 0)
+            .rotateX(direction.getTheta())
+            .rotateZ(direction.getPhi() + M_PI_2) *
         radius;
 
     // Get intersection with sphere
     const G4ThreeVector& toCenter = positionInDisk;
-    double t = sqrt(radius * radius - toCenter.dot(toCenter));
-    auto position = positionInDisk - t * direction;
+    const double t = sqrt(radius * radius - toCenter.dot(toCenter));
+    auto position = positionInDisk + t * direction;
 
     return position;
 }
@@ -237,8 +233,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event) {
             // radius in mm
             fCosmicCircumscribedSphereRadius = fSimulationManager->GetRestMetadata()
                                                    ->GetGeant4PrimaryGeneratorInfo()
-                                                   .GetSpatialGeneratorCosmicRadius() *
-                                               10.;
+                                                   .GetSpatialGeneratorCosmicRadius();
         }
 
         // This generator has correlated position / direction, so we need to use a different approach

src/SimulationManager.cxx

@@ -56,27 +56,29 @@ void PeriodicPrint(SimulationManager* simulationManager) {
                 outputPercentageType = "time";
             }
         }
-        const string eventsProgress =
-            " | " + to_string(simulationManager->GetNumberOfProcessedEvents()) + " events processed / " +
-            to_string(restG4Metadata->GetNumberOfEvents()) + " requested (" +
-            TString::Format(
-                "%.2e", simulationManager->GetNumberOfProcessedEvents() / simulationManager->GetElapsedTime())
-                .Data() +
-            "/s)";
-        const string timeProgress =
-            " | " + ToTimeStringLong(simulationManager->GetElapsedTime()) + " elapsed / " +
-            ToTimeStringLong(simulationManager->GetRestMetadata()->GetSimulationMaxTimeSeconds());
-
-        string progress = "";
+
+        string progress;
         if (outputPercentageType == "events") {
-            progress = eventsProgress;
+            progress = " | " + to_string(simulationManager->GetNumberOfProcessedEvents()) +
+                       " events processed / " + to_string(restG4Metadata->GetNumberOfEvents()) +
+                       " requested (" +
+                       TString::Format("%.2e", simulationManager->GetNumberOfProcessedEvents() /
+                                                   simulationManager->GetElapsedTime())
+                           .Data() +
+                       "/s)";
+        } else if (outputPercentageType == "entries") {
+            progress = " | " + to_string(simulationManager->GetNumberOfStoredEvents()) + " entries / " +
+                       to_string(restG4Metadata->GetNumberOfRequestedEntries()) + " requested";
         } else if (outputPercentageType == "time") {
-            progress = timeProgress;
+            progress = " | " + ToTimeStringLong(simulationManager->GetElapsedTime()) + " elapsed / " +
+                       ToTimeStringLong(simulationManager->GetRestMetadata()->GetSimulationMaxTimeSeconds());
         }
-        string timeInfo = "";
+
+        string timeInfo;
         if (outputPercentageType != "time") {
             timeInfo = " | " + ToTimeStringLong(simulationManager->GetElapsedTime()) + " elapsed";
         }
+
         G4cout << TString::Format("%5.2f", completionPercentage).Data() << "% | "
                << simulationManager->GetNumberOfStoredEvents() << " events stored / "
                << simulationManager->GetNumberOfProcessedEvents() << " processed ("