EvtGen does not take neutrino interaction probability into account

[sjoerd-bouma]

Currently, when a user specifies a particular neutrino energy spectrum in generator.py, the distribution of event energies matches this flux. In simulation.py, these events are then weighted by their survival probability after earth attenuation. However, their interaction probability is not taken into account. This means the resulting flux is not the one that would be observed by a detector, but instead matches the effective volume.

I think this is not particularly intuitive. I see two options:

1. We clearly document that this is the case, i.e. that the resulting event distribution after detector simulation matches the effective volume, and not the expected distribution at detection, and remind the user to apply some reweighting proportional to 1/ (interaction length) to obtain a spectrum that matches the expected distribution.
2. We fold the interaction probability with the specified neutrino energy spectrum in generator.py before sampling event energies (or at least add a parameter to enable this - I think this should be enabled by default).

My personal view on this is that 1. is error-prone, whereas 2. (if we enable this by default, which I think we should) changes the way we compute effective volumes or expected number of events in NuRadioMC.utilities.fluxes.py. The latter would still be my preference, but I'm open to discussion.

[fschlueter]

We have already started the discussion here: https://icecube-spno.slack.com/archives/G01GECH3GUX/p1707221347547539 (for context).

Okay, I do not have a strong opinion about this (yet). Let me just spill out a few thoughts:

• We should keep it simple. Right now it is quite simple, it happens in the step converting eff. volume to area. Adding a new weight factor at first does not sound terrible more complicated but it might becomes that once we add more weighting factors for what ever reason.
• The effective volume as it is defined now is independent of a model of the cross-section. Hence it is better when comparing detector performance. And it is somewhat a widely use observable in UHE neutrino telescopes. Changing that now causes at least confusion.
• We should avoid at all cost that it gets to easy to be confused an apply the energy dependent factor twice.

How would you do it @sjoerd-bouma, would you store a weight like L_int(E) / L_0 and when converting the effective volume you would include again the L_0?

[cg-laser]

Good that you bring this up. This is an important topic and we need to find a good solution. Probably best to have a call on that topic.

Here already a few thoughts:

• Felix raises good points that "interaction probability weights" will depend on the cross-section model
• I think the weight should be 1/Lint(E) so that Aeff = Veff * weight_interaction * weight_survival. (which is what we already do). We could just save such a weight in addition
• So far, we calculate event rate by using the Veff at specific energies (often 0.5 of 0.25 steps in logE). For each energy bin, we calculate the effective area and then the event rate by multiplying with an arbitrary flux, i.e. we approximate the integral via a sum. In case we would simulate a certain spectrum from the start, we would need to know what spectrum was simulated and then reweight the spectrum to the target spectrum. I think IceCube optical does it like this, so we should find out how their procedure is.
• The best approach will depend on the use case. For estimating detector sensitivity, it is best to simulate the Veff in fixed energy bins, and simulate enough events per bin to have a sufficiently large number of triggers. For generating DL training data sets, it is best to have an equal number of triggered events per energy bin, and then parameterize the uncertainty as a function of energy. For other approaches, it might be useful to have the event statistics represent their expected distribution, which I think is what Sjoerd wants to do.

[sjoerd-bouma]

So my suggestion in 2), if that wasn't clear, was sampling from f(E)/Lint(E) rather than sampling from f(E) as is done now. It's true that this introduces some (additional - of course, the CC/NC ratio and inelasticities already depend on the cross section) model dependence. But this would for example not change anything for simulations at fixed E.

[sjoerd-bouma]

To summarize the discussion in the inicemc call yesterday: we decided not to add another parameter named weight or probability due to the potential confusion with the pre-existing weight particle parameter which encodes the survival probability after attenuation in the Earth. Instead, the plan is to change the sampling in EvtGen to sample from f(E)/L_int(E) instead of f(E):

• generated events will follow ‘correct’ distribution automatically
• A_eff(E) = V_eff(E) / L_int(E) will still be valid at fixed E; doesn’t break compatibility with fluxes.py, which is only valid at fixed E anyway.

We will in addition store the interaction_length as a parameter, which will still allow for arbitrary reweighting of the flux afterwards.

[thoglu]

Hi, here is how IceCube does it. First there is a full Earth Monte Carlo being done to correctly simulate tau-regeneration or other absorption effects by the earth (using a specific Earth density model). Once those neutrinos reach the "active detection volume" (a region around IceCube which is sufficiently large to remove corner effects), we force an interaction. We always store the full individual weight information (interaction weight for forced interactions, zenith directional weight etc.) per event. But we also combine everything into a "oneweight" (which is also stored per event), which corresponds to the weight an E^-1 flux would have with a normalization of "1". Then you just need to multiply that "OneWeight" with a flux normalization an energy dependence if you do a reweighting for histograms. In principle you can also reconstruct the "oneweight" yourself. Typically, we genreate 50/50 anti neutrinos and neutrinos, but also that can be changed and corrected with appropriate "TypeWeights".

The whole procedure is described in the documentation:

Weights: https://software.icecube.wisc.edu/icetray/main/projects/neutrino-generator/weighting.html
Parameters in the WeightDict (Weight information per event): https://software.icecube.wisc.edu/icetray/main/projects/neutrino-generator/weightdict.html
How Neutrino generation works: https://software.icecube.wisc.edu/icetray/main/projects/neutrino-generator/intro.html

Finally, we have another project,
https://github.com/icecube/simweights,
which does the weighting for us, so we do not do an error while reading those per-event weight informations.. becaues it is quite easy to mess up an individual weight factor if you do the weighting yourself.
It also includes convenience functions to combine different datasets with different energy ranges etc.. this is the hard part,
because in order to calculate correct weights for combined datasets, you need to use all the individual weights (e.g. different datasets with different simulated zenith ranges, energy ranges and energy spectra), and combine them appropriately which requires some thinking.
I hope these ideas give some inspiration.

[fschlueter]

Maybe it is worth thinking about to also use simweights in NuRadio?

[thoglu]

yes I think that might be a good idea.. I remember 10 years ago in IceCube it was always a problem new people always need to understand the whole weighting scheme and everybody makes mistakes in the beginning, especially when combining datasets. It is very convenient to have a dedicated project that takes the weights as input and does the weighting for you and that everybody relies on.

[fschlueter]

Okay, I think for the short term I would support what @sjoerd-bouma summarised to be merged. Integrating simweights into NuRadio needs at least some carful testing... maybe more. I would volunteer looking into that if I find some time!

[thoglu]

I think the important part of the weight calculation is found here:
https://software.icecube.wisc.edu/icetray/main/projects/neutrino-generator/weighting.html#interaction-weight
The interaction is forced uniformly over the column depth in the active volume, and to correct for that and the forced interaction exponential weights (Interaction weight P_int - the probability to interact at all of the total column detph in active volume, pos-weight W_pos: we force the interaction uniformly, but in reality the interaction distribion is expontially distributed -see the table in the bottom) are saved that incorporate the exponential survival / interaction probabilities correctly. For each event, those values are differnet, but all of them are saved and the total weight can be recreated. Again, normally as a user one would just use the "OneWeight", as described above, and not worry about such details (or directly use the "simweights" tool).

[fschlueter]

I think the important part of the weight calculation is found here:
https://software.icecube.wisc.edu/icetray/main/projects/neutrino-generator/weighting.html#interaction-weight
The interaction is forced uniformly over the column depth in the active volume, and to correct for that exponential weights are saved that incorporate the exponential survival / interaction probabilities.

Good point! As far as I understand we do not take that into account (that the interaction probability varies over our detector volume)..

[sjoerd-bouma]

@thoglu Thanks for the input! I agree with Felix, having some tool similar to simweights for reweighting would be neat. In any case, it would be good to verify once more that we store everything that would be needed to compute a 'OneWeight', even if we don't want to explicitly do this (yet).

Re:

I think the important part of the weight calculation is found here:
https://software.icecube.wisc.edu/icetray/main/projects/neutrino-generator/weighting.html#interaction-weight
The interaction is forced uniformly over the column depth in the active volume, and to correct for that exponential weights are saved that incorporate the exponential survival / interaction probabilities.

Good point! As far as I understand we do not take that into account (that the interaction probability varies over our detector volume)..

- the interaction probability doesn't vary over the detector volume (other than the negligible but already included dependence on the survival probability). I'm pretty sure the interaction weight specified in the link is wrong (so I hope it's not implemented that way in the code @thoglu! Though I didn't check what effect it has). The survival probability P_survival is already the integrated survival probability, so to get the probability density of an interaction at X we need $\frac{d}{dX} (1-P_\mathrm{survival}) = \frac{1}{L_\mathrm{int}} P_\mathrm{survival}$ i.e. ignoring density fluctuations the only dependence on X is through the survival probability.

[fschlueter]

Okay I would agree that the interaction probability of an individual neutrino does not change over the volume. However, the spatial distribution (as a function of the direction) is not flat. But this is how we simulate it right?

[thoglu]

No @sjoerd-bouma , the derivation should be correct. There are two factors, so let us disentangle them.

First factor: probability of any interaction at all over the total length "L" in active volume.
The interaction probability (neglecting energy losses - which for neutrinos makes sense) per unit length is constant (Poisson), therefore the time/distance to interaction from the entry point has exponential form. The CDF over the exponential is the interaction probability, as we integrate the exponential. The CDF of an exponential is 1-exp(-C*x), and this is the weight. In the document they start already with the survival probability, which is 1-CDF, the "1s" cancel, and that is just the exponential factor. So in the document it is motivated the other way around (and might be confusing), but it still should be correct. That is the first factor, it is a probability (no density here!!!), and typically <<<<1.

The second factor is a relative weight of two "PDFs"! coming from the fact that we do not have uniform interaction distribution in reality (but of exponential form which is the exponential PDF of time/distance to interaction), but in simulation we "force" the interaction uniformly in IceCube, so we have a uniform distribution in simulation. So we have to correct for that by dividing the exponential PDF with a uniform PDF, a relative ratio of the two PDFs (which we call position weight). The second factor therefore is a ratio of two PDFs. The integral is there to properly normalize the exponential PDF.

The total weight is multiplication of both of those.

Where is the problem in the logic here? It is also in agreement with atmospheric and cosmic neutrino measurements by IceCube (in terms of the x,y,z distribution of measured neutrino events), which is further re-assuring.

[sjoerd-bouma]

@thoglu I agree with your explanation! But that's not what's in the text:

$$\begin{align}\begin{aligned}P_{pos\_True}(X) = 1 / sum * exp(- \sigma_{all} * X / M_{p} * C)\\\ sum = \int_0^{L_d} exp(- \sigma_{all} * X / M_{p} * C) dX\end{aligned}\end{align}$$

$P_{posTrue}$ is not the PDF. If the CDF is (1-P_survival(X)), the PDF is d/dX of this, which is what I stated in my previous comment. If we instead do the integral in sum, we get for the above:

$$\begin{equation} P_{posTrue} = P_{survival} / (- L_{int} [P_{survival}]^{L_D}_0) = \frac{P_{survival}}{L_{int} (1-P_{survival})} \end{equation}$$

Fortunately for IceCube, $P_{survival}\approx 1$, so we can expand the exponential and just get $1/L_D$ (i.e. the ratio of "PDFs" is just 1). Though I'm not sure where the dependence on $L_{int}$ has gone in that last result... so either I'm missing something or something is still wrong.

[thoglu]

@sjoerd-bouma P_PosTrue is a PDF. You and the document are both correct. The PDF the interaction probability is d/dX (P_int (X)=d/dX (1-exp(-Cx))) =Cexp(-C*X) which is just the exponential with a prefactor (that means 1-CDF=survival and PDF are related just by a prefactor, both of them have exponential form). However, in the document, the prefactor is not there (they only use the survival term for whatever reason as a "proportionality" to exp), so they calculate the prefactor numerically as 1/sum by forming the integral over the exponential. If you look at an exponential PDF, there has to be a prefactor (https://en.wikipedia.org/wiki/Exponential_distribution) for it to be correctly normalized. It is confusing because 1-CDF (survival) and PDF have both exponential form, but one has an extra prefactor.

I am pretty sure in the code there is no numerical integration, but just the respective analytical factors.

I think your analytic calculation is wrong though:
$\mathrm{SUM}=\int_{0}^{L} e^{-Q \cdot X} dX = \frac{1}{Q} \cdot (1- e^{-Q\cdot L}) \rightarrow P_{pos_True}(X)=\frac{exp(-Q \cdot X)}{\frac{1}{Q} \cdot(1- e^{-Q \cdot L})} \approx exp(-Q \cdot X)/L$ for small $Q \cdot L$ and further $\approx \frac{1-Q \cdot X}{L} \approx 1/L$ for small $Q \cdot X$, which however is already implied by $Q \cdot L$ small , where $Q=\sigma_{all}\cdot C / M_p$

so the ratio to 1/L, which is the uniform PDF, is almost 1 but not quite and that is the second weight

[sjoerd-bouma]

@thoglu Okay, I think I'm convinced now : ) took me long enough, but it seems they are indeed equivalent. Thanks for your patience!

I think my analytic calculation is right and yours is wrong though ;) the 1/Q goes outside. It has units, so your expression doesn't work

[thoglu]

corrected, but I don't see them to be the same still

[sjoerd-bouma]

Yeah, I was a bit sloppy with my notation - the $P_{survival}$ on top should be $exp(-QX)$ and the one on the bottom $exp(-QL)$. You are of course correct that this is a (very small!) correction to the first factor. But yes, in principle it is necessary. My confusion arose mostly because we keep track of this term differently - we define the survival probability directly as $exp(-QX)$ but with $X$ the propagated distance through the entire Earth until the interaction point (i.e. we don't separately calculate the survival probability until the edge of the detector volume or the interaction probability within the detector volume), so the dependence on the interaction length shows up differently in both formulas.

$$\begin{equation} P_{posTrue} = P_{survival}(X) / (- L_{int} [P_{survival}]^{L_D}_0) = \frac{P_{survival}(X)}{L_{int} (1-P_{survival}(L_D))} \end{equation}$$

The first factor was just $P_{interaction} = 1- P_{survival}(L_D)$ so multiplying these together just gives you the same result I gave in my earlier comment. (Modulo my sloppy notation/misunderstanding of what $P_{survival}$ means in your documentation vs how it is defined in NuRadioMC)

[thoglu]

yeah in IceCube neutrinos are propagated through the Earth, and if they do not make it to the active volume, they are dropped.. there is no forced interaction for neutrinos that do not reach the active volume (in default mode at least)

[cozzyd]

I don't know if this is useful at all or not, but I think this related to a point I brought up several years ago? Basically the uniform volume sampling is not really correct for calculating A_eff, I think.
Aeff_slides.pdf
(in these slides, I'm ignoring survival weight, just assume that you throw out neutrinos that don't make it there, but you could equivalently multiply p_obs by survival weight).

[thoglu]

It depends I guess how you use it. In IceCube, we never really work with V_eff, for example. We directly generate a flux from a "generating surface", not a "generating volume", and track everything with respect to that surface to translate to fluxes etc. the uniform sampling over the length in the volume then just corrects with a given weight and should be correct.

[cozzyd]

Yeah what IceCube does (work from a generated surface) seems more sensible, since V_eff can't be exactly translated to A_eff except in a few special geometries, but we could store the per-event conversion factor as some sort of interaction weight (that's slightly more complicated than 1/l_int; I think it's more like A_proj(theta)/( * l_int or something like that, as long as l_int doesn't get too small).

[cg-laser]

Hi Cosmin, all. I only vaguely remember this discussion from three years ago but back then I went through the effort of implementing an unforced event generator to verify that we indeed simulate things correctly with our "forced interaction" method. See slides here:
20200325_forced_unforced_evtgen.pdf

[cozzyd]

Yes I remember this now. Can you remind me if you're just counting events within cylinder (i.e. the trigger condition is that the vertex is inside the cylinder, or something similar) or if there is a more complicated "trigger" condition? I think for the methods to diverge, the triggering must depend on the spatial position and neutrino direction in a non-uniform way (for simple trigger conditions you will get the same answer since the position and direction can be factored apart).

I tried something similar (where I had a simple emission model + straight line raytracing) and got a ~7 percent difference in diffuse effective areas at high energies (See https://github.com/cozzyd/ToyDisc/blob/master/slides/ToyDisc.pdf ). So probably not a big effect, but it seemed to be there. There could be a bug of course. But by adding some sort of interaction weight, I think it could be verified with a more realistic simulation?

[sjoerd-bouma]

@cozzyd I've been wrong before in this thread, but I don't think there's anything wrong with the current simulation method.

The interaction PDF is still $\frac{d}{dX} (1-P_\mathrm{survival}) = \frac{1}{L_\mathrm{int}} P_\mathrm{survival}$. Imagine splitting up the cylinder (or indeed an arbitrary simulation volume) into tiny cubic pixels (voxels?) of side $l$. Then for an isotropic flux $I_0$, the probability of an interaction inside this pixel is

$$I \times \int_{A_{pix}} \frac{l_{slant}}{L_{int}} \times P_\mathrm{survival} \cdot dA = I \times \int_{A_{pix}} (l_{slant} \cdot dA) \times \frac{P_\mathrm{survival}}{L_{int}}$$

(if we're interested in the number of triggering interactions, we just multiply by $p_{obs}$, but this 'weighting' happens automatically in the simulation). But $\int_{A_{pix}} l_{slant} \cdot dA$ is just, as you state in different words in your slides, the volume of the pixel, because a volume doesn't depend on how you orient your axes when you integrate it. So each pixel just contributes $V_{pix} / L_{int} \cdot P_\mathrm{survival}$, i.e. the number of interactions is just proportional to the volume of each pixel, not its projected area (well, technically through $L_{int}$ it's proportional to the number of nucleons in the pixel). In particular, the only dependence on $\theta$ is through $P_\mathrm{survival}$ (because again, $p_{obs}$ is included automatically simply by running the trigger simulation).

The long and short of it is as long as we take into account the factors $1/L_{int} \times P_\mathrm{survival}$ we are correct to simulate uniformly in the volume.

[cozzyd]

So I wonder if the difference is entirely due to the exponential sampling of the interaction point along the trajectory when I use the area method. That would explain why it's a bigger effect at higher energies, and it also breaks the "each pixel is independent of each other" argument, while remaining a smallish effect in the in-ice case (while being a big effect in the balloon case where we first noticed the difference!). If this is the case, then it's easy to fix with an interaction weight exactly as Thorsten described (the ratio of the exponential to uniform PDF along the intersecting segment). I can test this pretty easily with my Toy MC by sampling linearly rather than exponentially and seeing if the two sampling methods then agree.

(Ok, this would still neglect that the ice density is not uniform, but that can be incorporated too).

[sjoerd-bouma]

So I wonder if the difference is entirely due to the exponential sampling of the interaction point along the trajectory when I use the area method. That would explain why it's a bigger effect at higher energies, and it also breaks the "each pixel is independent of each other" argument, while remaining a smallish effect in the in-ice case (while being a big effect in the balloon case where we first noticed the difference!). If this is the case, then it's easy to fix with an interaction weight exactly as Thorsten described (the ratio of the exponential to uniform PDF along the intersecting segment). I can test this pretty easily with my Toy MC by sampling linearly rather than exponentially and seeing if the two sampling methods then agree.

(Ok, this would still neglect that the ice density is not uniform, but that can be incorporated too).

So maybe I'm missing something, but the argument I gave in my previous comment is not approximate, it is exact. One thing to note is that our definition of $P_\mathrm{survival}$ is different from the one used by IceCube, i.e. it is the survival probability to the neutrino vertex and doesn't distinguish between attenuation inside and outside the simulation volume. What we are probably doing wrong is just using PREM for the survival probability rather than forcing the outside layers to be ice and bedrock, but (I haven't checked) PREM probably assumes something like the density of water, which is not too far off.

[sjoerd-bouma]

There really aren't separate survival probabilities inside and outside the simulation volume, and neither is there an 'interaction probability' separate from the survival probability!

$$\frac{d}{dx}|_{\vec{x}=\vec{x}_\mathrm{int}} (1-P_\mathrm{survival}(\vec{x})) = \frac{1}{L_\mathrm{int}(\vec{x}_\mathrm{int})} P_\mathrm{survival}(\vec{x}_\mathrm{int})$$

for the interaction PDF is completely general, even the dependence on the ice density is encoded in $L_\mathrm{int}$

[cozzyd]

So I'm not sure I agree that you can specify that you split up the volume into voxels without specifying their shape, I think you may just end up with the same problem on a smaller scale?

Based on numerical experiments, the exponential placement of the vertex indeed has nothing to do with anything (it's a good enough approximation, at least for the size volume considered here).

The problem I think that for most geometries (including a cylinder), the projected area is different from different angles, which (I think) means that it's not appropriate to sample uniformly in direction if sampling volumetric. For example, here is the distribution of direction Z component as a function of distance^2 from center of the cylinder with the volumetric and area-based sampling methods:

The area-based sampling method has to be weighted by interaction probability (otherwise the corners "catch" a lot more events, but since the path length is shorter, you're less likely to interact there; this is using a 1 EeV cross-section. The reason I show the perp2 component is so you can see there is no incorrect pileup in the corners with the area-based method).

In both cases, the distribution of vertices is indeed uniform within the volume, but the angular distribution of arriving neutrinos is not, even before applying the survival probability. There could be a bug (or many!) in my toy MC but I think the idea makes intuitive sense.

As for PREM, yes, you probably want to modify PREM so that the top 3 km (or whatever) are ice. It kind of looks like that's almost the case, except NuRadioMC using 1.002 g/cm^3 instead of 0.917 g/cm^3 for the top layer, so it's probably not too important. (My toy MC has it match, otherwise there would be strange discontinuities related to the fact that the survival probability for the area-based sampling is defined at entrance into the volume). There will still be a discontinuity in the firn using a more realistic ice model, but that's probably even more negligible...

Though if a significant fraction of events do happen in the firn (which I think has to be true for some energies /zenith angles), then using l_int based on deep ice is also not going to be right. Since the density of firn varies by more than a factor of 2, that might actually be pretty important to get right...

[sjoerd-bouma]

I did specify the shape of my voxels! They're cubic with edge length $l$. But you just end up integrating the projected area times the slant depth, which as I've stated is just the volume, so the shape doesn't matter as long as you make them small enough (because then the survival probability and the detection probability factorize).

It's not the projected area which is important, but the volume. This makes sense! In the case where the survival probability can be ignored, it's clear (I hope) that all that matters is the number of nucleons you have, and in fact this approximation is probably already better than the 7% discrepancy you found in your toy MC. If the survival probability can not be ignored, the result is exactly the one I stated previously, which is how it is treated already (okay, except PREM not having the correct density at the edges and as you mention $L_{int}$ should be position-dependent as well).

[sjoerd-bouma]

The area-based sampling method has to be weighted by interaction probability (otherwise the corners "catch" a lot more events, but since the path length is shorter, you're less likely to interact there; this is using a 1 EeV cross-section. The reason I show the perp2 component is so you can see there is no incorrect pileup in the corners with the area-based method).

This 'area vs path length' trade-off is exactly why you just get $\int l_{slant} \cdot dA$ which is just the volume. It is easier to see when you split the volume up into voxels because then $p_{obs}$, $L_{int}$ and $P_\mathrm{survival}$ can just be taken outside the integral. So as long as you take all of them into account after sampling randomly inside your volume you will get the correct distribution.

[sjoerd-bouma]

Another reason why the projected area doesn't matter - in the extreme case you mention ('perfect detection efficiency'), all/most neutrinos would interact essentially as soon as they hit the detector volume (or most likely before), so that the survival probability for anything that has travelled a significant distance within the detector volume goes to 0. In that case, uniform volume sampling * survival probability gets you... something that looks very much like the projected area, as you would expect!

[cozzyd]

So I agree that the number of interactions just depends on number of nucleons, but in our case, we also care about the angular distribution of interacting particles, which will affect what you can observe (with an Askaryan detector, at least... with a bucket of scintillator, that's not something that affects what you observe). So I guess it's wrong to say that there's a problem with the volumetric sampling of vertex position, but instead maybe there's a problem with the isotropic sampling of neutrino direction within the volume, since the volume will have preferred directions (even if it was put out in space and there was no survival probability to talk about) due to having different cross-sections in different directions?

I certainly agree that the integral of
$l_{slant} \cdot dA$ is the volume, but it doesn't follow to me that you can necessarily separate it out from the rest of the integral?

Note also that the "perfect detection efficiency" doesn't mean that every neutrino interacts, but rather than you observe every neutrino that interacts. In that case, since p_{obs} does not depend on angle, maybe you can indeed integrate out the volume and get the same answer.

But I will be the first to admit that my geometric intuition regarding these things is not great, and I could be talking nonsense. (This is why I tried to do this toy MC thing... but there are lots of ways I could have screwed up...).

[sjoerd-bouma]

I certainly agree that the integral of lslant⋅dA is the volume, but it doesn't follow to me that you can necessarily separate it out from the rest of the integral?

For a small enough voxel, you can, because $p_{obs}$, $L_{int}$ and $P_\mathrm{survival}$ don't vary over the volume of the voxel. But if you do this each voxel contributes as $V$ times these factors, i.e. we can just sample uniformly over the entire volume as long as we include $p_{obs}$, $L_{int}$ and $P_\mathrm{survival}$ for each sampled vertex position (or voxel). There is still some dependence on the angle, but this is fully encoded into $P_\mathrm{survival}$.