docs: def scenario_damage(riskinputs, riskmodel, param, monitor):
      """
      Core function for a damage computation.
  
      :param riskinputs:
          :class:`openquake.risklib.riskinput.RiskInput` objects
      :param riskmodel:
          a :class:`openquake.risklib.riskinput.CompositeRiskModel` instance
      :param monitor:
          :class:`openquake.baselib.performance.Monitor` instance
      :param param:
          dictionary of extra parameters
      :returns:
          a dictionary {'d_asset': [(l, r, a, mean-stddev), ...],
                        'd_event': damage array of shape R, L, E, D,
                        'c_asset': [(l, r, a, mean-stddev), ...],
                        'c_event': damage array of shape R, L, E}
  
      `d_asset` and `d_tag` are related to the damage distributions
      whereas `c_asset` and `c_tag` are the consequence distributions.
      If there is no consequence model `c_asset` is an empty list and
      `c_tag` is a zero-valued array.
      """
      L = len(riskmodel.loss_types)
      D = len(riskmodel.damage_states)
      E = param['number_of_ground_motion_fields']
      R = riskinputs[0].hazard_getter.num_rlzs
      result = dict(d_asset=[], d_event=numpy.zeros((E, R, L, D), F64),
                    c_asset=[], c_event=numpy.zeros((E, R, L), F64))
      for ri in riskinputs:
          for out in riskmodel.gen_outputs(ri, monitor):
              r = out.rlzi
              for l, loss_type in enumerate(riskmodel.loss_types):
                  for asset, fractions in zip(ri.assets, out[loss_type]):
                      dmg = fractions[:, :D] * asset['number']  # shape (E, D)
                      result['d_event'][:, r, l] += dmg
                      result['d_asset'].append(
                          (l, r, asset['ordinal'], scientific.mean_std(dmg)))
                      if riskmodel.consequences:
                          csq = fractions[:, D] * asset['value-' + loss_type]
                          result['c_asset'].append(
                              (l, r, asset['ordinal'], scientific.mean_std(csq)))
                          result['c_event'][:, r, l] += csq
      return result