From 30eb90099b1fc5b6509cb4b202c8fa721c5971cd Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Tue, 14 Sep 2021 12:16:53 +0200
Subject: [PATCH] Introducing cn as parameter for the Expiration, and proper
 normalization of the exposure

---
 cara/models.py           | 15 ++++++++++++---
 cara/monte_carlo/data.py | 22 +++++++---------------
 2 files changed, 19 insertions(+), 18 deletions(-)

diff --git a/cara/models.py b/cara/models.py
index dd797ba2..18c3990a 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -922,9 +922,18 @@ class ExposureModel:
         return normed_exposure * self.repeats
 
     def exposure(self) -> _VectorisedFloat:
-        """The number of virus per meter^3."""
-        return (np.array(self._normed_exposure()).mean() *
-                self.concentration_model.infected.emission_rate_when_present())
+        """
+        The number of virus per meter^3. With sampled diameters, the
+        aerosol volume has to be put back in the exposure before taking
+        the mean, to obtain the proper result for the exposure (which
+        corresponds to an integration on diameters).
+        """
+        mask = self.concentration_model.infected.mask
+        aerosols = self.concentration_model.infected.expiration.aerosols(mask)
+        emission_rate = self.concentration_model.infected.emission_rate_when_present()
+
+        return (np.array(self._normed_exposure()*aerosols).mean() *
+                emission_rate/aerosols)
 
     def infection_probability(self) -> _VectorisedFloat:
         exposure = self.exposure()
diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py
index 09dd50d2..e033a558 100644
--- a/cara/monte_carlo/data.py
+++ b/cara/monte_carlo/data.py
@@ -52,14 +52,6 @@ class BLOmodel:
                     for A,cn,mu,sigma in zip(self.BLO_factors, self.cn,
                                              self.mu, self.sigma) )
 
-    def normalized_distribution(self, d, dmin, dmax):
-        """
-        Return the probability distribution, normalized by its integral
-        between dmin and dmax (microns).
-        """
-        norm = self.integrate(dmin, dmax)
-        return self.distribution(d) / norm
-
     def integrate(self, dmin, dmax):
         """ 
         Returns the integral between dmin and dmax (in microns) of the 
@@ -143,17 +135,17 @@ def expiration_distribution(BLO_factors: typing.Tuple[float, float, float]):
     """
     Returns an Expiration with an aerosol diameter distribution, defined
     by the BLO factors.
-    Note: integration boundaries for normalization are chosen as 0.1 and
-    30 microns respectively - this is an historical choice based on
-    previous implementations of the model (it limits the influence of
-    the O-mode).
+    The total concentration of aerosols is computed by integrating
+    the distribution between 0.1 and 30 microns - these boundaries are
+    an historical choice based on previous implementations of the model
+    (it limits the influence of the O-mode).
     """
+    dscan = np.linspace(0.1, 30. ,3000)
     return mc.Expiration(CustomKernel(dscan,
-                BLOmodel(BLO_factors).normalized_distribution(dscan, 0.1, 30.),
-                kernel_bandwidth=0.1))
+                BLOmodel(BLO_factors).distribution(dscan),kernel_bandwidth=0.1),
+                BLOmodel(BLO_factors).integrate(0.1, 30.))
 
 
-dscan = np.linspace(0.1, 30. ,3000)
 expiration_distributions = {
     'Breathing': expiration_distribution((1., 0., 0.)),
     'Talking':   expiration_distribution((1., 1., 1.)),