Closed reichherzerp closed 2 years ago
Thanks for the hint, hmm ok! Great that you could improve the code!
I just inherited the TabularPhotonField class to the PhotonFieldSampling class to use its functionality there and call the constructor of TabularPhotonField therefore now in every instanciation of the PhotonFieldSampling constructor. This way the tabFiles are read and the photonEnergies array is directly available in sample_eps. However, I was not sure when calling the constructor for the inheritance with a certain parameter set how to call it afterwards again selecting the correct field depending on the bgFlag. Therefore, I only solved it a bit unelegantly and in my first attempt now added the flag selecting the photonfield to the TabularPhotonField constructor (I think the TabularPhotonField class should not use the flag and is already more elaborated using fieldnames). However, I think it works and I checked that it should read the correct energies.
I think I found out why the file containing the densities is so much longer than the file containing the energies, it contains the densities for the different redshifts given in the redshift file. length(tabDensity) = length(tabEnergy) * length(tabRedshift)
Thanks for the improvements! It is great that It works as you implemented it so that we can finalize the other parts of the sampling already. We can discuss the tabFiles reading process on Monday. For now, I will continue working on the improved calculation of the pMax by using interpolated photon energies from the tab files. I'll also remove the computation of the mass and the momentum as we don't need them anymore for the sampling.
We could think about warning the user if the provided custom fields have a "wrong" energy range for which no photon satisfies the criteria for interaction. In this case, we would be stuck in the sampling loop forever:
This is important, as CRPropa has decided already that we have to get a photon when we enter the sampling function! When we fail to do so, we have to warn the user...
After thinking some more time about it, I think that we could implement 2. as a hard criterion. Here, we can throw dedicated and detailed error messages, so that the user knows immediately what to do. However, there may be even cases where the 2. criteria is just met, but the sampling still would need too many iterations. For this, I vote for having the 1. option with a more generic error message as well. Looking forward to your thoughts on this.
@JulienDoerner, @LeanderSchlegel: I noticed that we have to calculate the maximum probability pMax in every interaction because it depends on the current particle energy. We can't compute it only once in the beginning...
However, I noticed that we do not have to calculate cnorm at all (both values in the MC comparison are divided by cnorm, which is the same as to only compare the initial values) and since I have removed it now the code should already become faster!