Integration of solar axion flux

[jgalan]

• TRestAxionSolarModel is kept as an interface to sebhoof/SolarAxionFlux repository.
• TRestAxionSolarFlux has been recovered from the previous TRestAxionSolarModel implementation. Now it is possible to load a continuum spectra component, and a monochromatic component.
• Implemented method TRestAxionSolarFlux::GetRandomEnergyAndRadius that returns a random energy and position following the distribution from the tables loaded in the object.
• TRestAxionSolarFlux has been fully documented.
• A script has been produced to generate few plots and test the TRestAxionSolarFlux class.

We need yet to implement it at TRestAxionGeneratorProcess to produce random Monte Carlo events and process them in a full event chain.

[jgalan]

The following plots are generated with the script that is now executed at the pipeline.

[jgalan]

Hi there, TRestAxionSolarFlux has been fully documented, there is a new script that can be used as an example to produce the plot presented in this PR, and it has been integrated in the pipeline for validation.

The corresponding axion generator process is also available at #12

@sebhoof @jovoy @rest-for-physics/axionlib if there you dont's see a major or minor issue, plz, approve this PR

[sebhoof]

The PR adds the desired features. The TRestAxionSolarFlux::GetRandomEnergyAndRadius() function draws random radii and energy values from the 2D distribution of continuum and monochromatic lines.

I did not find any major issues (however, see comment below), so the PR can be merged IMO.

Comment for the future: the implementation effectively assumes a step function interpolation of the underlying radial distribution (since the tabulated data contains the integrated flux in rings on the solar disc). The rings thus need to be small (the current value of dr = 0.01 Rsol is fine) because the underlying r-distribution can never be recovered as a continuous function. Another reason is that artificially many axions are generated with r ~ 0 in the 1st bin only; again, this is not an issue as long as dr is small.

[jgalan]

Hi @sebhoof we will not be able to merge if there is no positive review (at least one positive review is required). I added you to the reviewers as the right candidate to approve merging of this PR. But the merging will be blocked till you do not explicitly approve it.

[jgalan]

Once you approve you will be able to merge yourself, I believe.

[jgalan]

Another reason is that artificially many axions are generated with r ~ 0 in the 1st bin only; again, this is not an issue as long as dr is small.

Yes, it is not the case, as it can be seen on the figure showing the MonteCarlo counting versus the solar radius. If the tables are properly integrated I dont see how there could be an artificial excess, but of course, there would be a lost on spatial solar resolution.

[sebhoof]

Yes, sorry, I forgot to click the submit review button after doing some more tests.

The artificial excess can happen since no axions with r = 0 should be generated; you're right that this can be understood as a loss in spatial resolution as well. Instead of assigning an event in radial bin i to a fixed radius, you could draw another random number in U(0,1) and map the event to some random radius within radial bin i.

In any case, as long as the radial bins are small enough, this is not a problem. It's just important to note that "small enough" does not depend on the actual spatial resolution of the optics/detector but due to the loss of information introduced in the MC scheme.

[jgalan]

The artificial excess can happen since no axions with r = 0 should be generated; you're right that this can be understood as a loss in spatial resolution as well. Instead of assigning an event in radial bin i to a fixed radius, you could draw another random number in U(0,1) and map the event to some random radius within radial bin i.

Yes, I believe what you describe is what we have at GetRandomEnergyAndRadius.

                Double_t energy = fFluxTable[r]->GetRandom();
                 Double_t radius = ((Double_t)r + fRandom->Rndm()) * 0.01;

fFluxTable is the spectrum at a given radius, the energy returned by TH1D object I am sure is not discretized with the binning but each bin energy is mapped using U(0,1) uniform distribution.