Fixing solar axion generator

.github/pr-badge.yml

@@ -17,3 +17,5 @@
   url: "https://gitlab.cern.ch/rest-for-physics/axionlib/-/commits/$branchName"
 - imageUrl: "https://gitlab.cern.ch/rest-for-physics/framework/badges/$branchName/pipeline.svg"
   url: "https://gitlab.cern.ch/rest-for-physics/framework/-/commits/$branchName"
+- imageUrl: "https://github.com/rest-for-physics/axionlib/actions/workflows/validation.yml/badge.svg?branch=$branchName"
+  url: "https://github.com/rest-for-physics/axionlib/commits/$branchName"

inc/TRestAxionGeneratorProcess.h

@@ -57,11 +57,8 @@ class TRestAxionGeneratorProcess : public TRestEventProcess {
     /// Seed used in random generator
     Int_t fSeed = 0;  //<
 
-    /// It defines the minimum energy as a cut-off to the generator. Default is 50eV.
-    Double_t fMinEnergy = 0.05;  //<
-
-    /// It defines the maximum energy as a cut-off to the generator. No limits if equals 0.
-    Double_t fMaxEnergy = 0;  //<
+    /// It defines the energy range for the axion event generator. Default between 50eV and 10keV.
+    TVector2 fEnergyRange = TVector2(0.05, 10);  //<
 
     void Initialize() override;
 

inc/TRestAxionSolarFlux.h

@@ -35,19 +35,17 @@
 //! A metadata class to load tabulated solar axion fluxes
 class TRestAxionSolarFlux : public TRestMetadata {
    private:
-    void Initialize();
-
     /// The filename containning the solar flux table with continuum spectrum
     std::string fFluxDataFile = "";  //<
 
     /// The filename containning the solar flux spectra for monochromatic spectrum
     std::string fFluxSptFile = "";  //<
 
     /// Axion coupling. Defines coupling type and strength.
-    std::string fCouplingType;  //<
+    std::string fCouplingType = "";  //<
 
     /// Axion coupling strength
-    Double_t fCouplingStrength;  //<
+    Double_t fCouplingStrength = 0;  //<
 
     /// Seed used in random generator
     Int_t fSeed = 0;  //<
@@ -97,6 +95,8 @@ class TRestAxionSolarFlux : public TRestMetadata {
     /// A metadata member to control if the tables have been loaded
     Bool_t fTablesLoaded = false;  //!
 
+    void Initialize();
+
     void ReadFluxFile();
     void LoadContinuumFluxTable();
     void LoadMonoChromaticFluxTable();
@@ -112,7 +112,7 @@ class TRestAxionSolarFlux : public TRestMetadata {
     /// It returns the total integrated flux at earth in cm-2 s-1
     Double_t GetTotalFlux() { return fTotalContinuumFlux + fTotalMonochromaticFlux; }
 
-    std::pair<Double_t, Double_t> GetRandomEnergyAndRadius();
+    std::pair<Double_t, Double_t> GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1));
 
     void LoadTables();
 

src/TRestAxionGeneratorProcess.cxx

@@ -177,11 +177,10 @@ TRestEvent* TRestAxionGeneratorProcess::ProcessEvent(TRestEvent* evInput) {
 
     r = TMath::Sqrt(r);
 
+    std::pair<Double_t, Double_t> p;
     if (fGeneratorType == "solarFlux") {
-        std::pair<Double_t, Double_t> p = fAxionFlux->GetRandomEnergyAndRadius();
+        p = fAxionFlux->GetRandomEnergyAndRadius(fEnergyRange);
         energy = p.first;
-        if (energy < fMinEnergy) return nullptr;
-        if (energy > fMaxEnergy) return nullptr;
         Double_t radius = p.second;
 
         axionPosition = TVector3(REST_Physics::solarRadius * radius * x / r,
@@ -191,10 +190,13 @@ TRestEvent* TRestAxionGeneratorProcess::ProcessEvent(TRestEvent* evInput) {
     }
 
     if (fGeneratorType == "flat" || fGeneratorType == "plain") {
-        if (fMaxEnergy > 0)
-            energy = fRandom->Rndm() * (fMaxEnergy - fMinEnergy) + fMinEnergy;
-        else
-            energy = (1. - fMinEnergy) * fRandom->Rndm() + fMinEnergy;
+        if (fEnergyRange.X() > 0 && fEnergyRange.Y() > 0)
+            energy = fRandom->Rndm() * (fEnergyRange.Y() - fEnergyRange.X()) + fEnergyRange.X();
+        else if (fEnergyRange.X() > 0)
+            energy = (1. - fEnergyRange.X()) * fRandom->Rndm() + fEnergyRange.X();
+        else {
+            RESTWarning << "Not a valid energy range was defined!" << RESTendl;
+        }
     }
 
     /// The axion position must be displaced by the target size.
@@ -232,7 +234,7 @@ void TRestAxionGeneratorProcess::PrintMetadata() {
     RESTMetadata << "Axion mass: " << fAxionMass * units("eV") << " eV" << RESTendl;
     RESTMetadata << "Target radius: " << fTargetRadius * units("cm") << " cm" << RESTendl;
     RESTMetadata << "Random seed: " << (UInt_t)fSeed << RESTendl;
-    RESTMetadata << "Energy range: (" << fMinEnergy << ", " << fMaxEnergy << ") keV" << RESTendl;
+    RESTMetadata << "Energy range: (" << fEnergyRange.X() << ", " << fEnergyRange.Y() << ") keV" << RESTendl;
 
     RESTMetadata << "+++++++++++++++++++++++++++++++++++++++++++++++++" << RESTendl;
 }

src/TRestAxionSolarFlux.cxx

@@ -682,7 +682,7 @@ void TRestAxionSolarFlux::IntegrateSolarFluxes() {
 /// \brief It returns a random solar radius position and energy according to the
 /// flux distributions defined inside the solar tables loaded in the class
 ///
-std::pair<Double_t, Double_t> TRestAxionSolarFlux::GetRandomEnergyAndRadius() {
+std::pair<Double_t, Double_t> TRestAxionSolarFlux::GetRandomEnergyAndRadius(TVector2 eRange) {
     std::pair<Double_t, Double_t> result = {0, 0};
     if (!fTablesLoaded) return result;
     Double_t rnd = fRandom->Rndm();
@@ -691,6 +691,9 @@ std::pair<Double_t, Double_t> TRestAxionSolarFlux::GetRandomEnergyAndRadius() {
         for (int r = 0; r < fFluxTableIntegrals.size(); r++) {
             if (rnd < fFluxTableIntegrals[r]) {
                 Double_t energy = fFluxTable[r]->GetRandom();
+                if (eRange.X() != -1 && eRange.Y() != -1) {
+                    if (energy < eRange.X() || energy > eRange.Y()) return GetRandomEnergyAndRadius(eRange);
+                }
                 Double_t radius = ((Double_t)r + fRandom->Rndm()) * 0.01;
                 std::pair<Double_t, Double_t> p = {energy, radius};
                 return p;