DDDigi/plugins/DigiDepositSmearEnergy.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include <DDDigi/DigiContainerProcessor.h>
0016 #include <DD4hep/DD4hepUnits.h>
0017
0018 /// C/C++ include files
0019 #include <limits>
0020
0021 /// Namespace for the AIDA detector description toolkit
0022 namespace dd4hep {
0023
0024   /// Namespace for the Digitization part of the AIDA detector description toolkit
0025   namespace digi {
0026
0027     /// Actor to smear geographically energy deposits
0028     /**
0029      *   Intrinsic fluctuation:     sigma(E)/E = const(fluctuation) / sqrt(E)
0030      *                              sigma(E)   = const(fluctuation) * sqrt(E)
0031      *   Sampling fluctuation:      sigma(E)/E = const(sampling) * sqrt(deltaE/E)
0032      *                              sigma(E)   = const(sampling) * sqrt(deltaE*E)
0033      *                              deltaE = energy loss to create single electron/ion
0034      *                                       pair in sampling layer
0035      *   Instrumetation:            sigma(E)/E = const(instrumentaion) / E
0036      *                              sigma(E)   = const(instrumentaion)
0037      *   Systematics:               sigma(E)/E = const(systematics)
0038      *                              sigma(E)   = const(systematics) * E
0039      *
0040      *   electron conversion
0041      *   statictical fluctuation:   poissonian change of converted e-ion pairs
0042      *
0043      *   Electromagnetic showers:
0044      *   Shower counters:
0045      *   Fluctuation resolution: sigma(E)/E ~ 0.005 / sqrt(E)
0046      *
0047      *   Sampling calorimeters:
0048      *   Sampling resolution:    sigma(E)/E ~ 0.04 * sqrt(1000*dE/E)
0049      *   dE = energy loss of a single charged particle in one sampling layer
0050      *
0051      *   Hadronic showers:
0052      *   Shower counters:
0053      *   Fluctuation resolution: sigma(E)/E ~ 0.45 / sqrt(E)
0054      *   With compensation for nuclear effects: sigma(E)/E ~ 0.25 / sqrt(E)
0055      *   Sampling resolution:    sigma(E)/E ~ 0.09 * sqrt(1000*dE/E)
0056      *   dE = energy loss of a single charged particle in one sampling layer
0057      *
0058      *   (See all: https://handwiki.org/wiki/Physics:Energy_resolution_in_calorimeters)
0059      *
0060      *   Properties:
0061      *   intrinsic_fluctuation:      Intrinsic energy resolution constant (gaussian ~ std::sqrt(energy))
0062      *   systematic_resolution:      Systematic energy resolution constant (gaussian ~deposit energy)
0063      *   instrumentation_resolution: Instrumentation energy resolution constant (gaussian(CONSTANT))
0064      *   pair_ionization_energy:     Ionization constant to create electron ION pair in absorber
0065      *   ionization_fluctuation:     Flag to use ionization fluctuation smearing (poisson ~ # e-ion-pairs)
0066      *   modify_energy:              Flag to modify energy deposit according to error - otherwise only compute error
0067      *
0068      *  \author  M.Frank
0069      *  \version 1.0
0070      *  \ingroup DD4HEP_DIGITIZATION
0071      */
0072     class DigiDepositSmearEnergy : public DigiDepositsProcessor   {
0073     protected:
0074       /// Property: Intrinsic energy resolution constant (gaussian ~ std::sqrt(energy))
0075       double m_intrinsic_fluctuation      { 0e0 };
0076       /// Property: Systematic energy resolution constant (gaussian ~deposit energy)
0077       double m_systematic_resolution      { 0e0 };
0078       /// Property: Instrumentation energy resolution constant (gaussian(CONSTANT))
0079       double m_instrumentation_resolution { 0e0 };
0080       /// Property: Ionization constant to create electron ION pair in absorber
0081       double m_pair_ionization_energy     { 0e0 };
0082       /// Property: Flag to use ionization fluctuation smearing (poisson ~ # e-ion-pairs)
0083       bool   m_ionization_fluctuation     { false };
0084       /// Property: Flag to modify energy deposit according to error - otherwise only compute error
0085       bool   m_modify_energy              { true };
0086
0087     public:
0088       /// Standard constructor
0089       DigiDepositSmearEnergy(const DigiKernel& krnl, const std::string& nam)
0090         : DigiDepositsProcessor(krnl, nam)
0091       {
0092         declareProperty("intrinsic_fluctuation",      m_intrinsic_fluctuation = 0e0);
0093         declareProperty("instrumentation_resolution", m_instrumentation_resolution = 0e0);
0094         declareProperty("systematic_resolution",      m_systematic_resolution  = 0e0);
0095         declareProperty("pair_ionisation_energy",     m_pair_ionization_energy = 3.6*dd4hep::eV);
0096         declareProperty("ionization_fluctuation",     m_ionization_fluctuation = false);
0097         declareProperty("modify_energy",              m_modify_energy = true);
0098         DEPOSIT_PROCESSOR_BIND_HANDLERS(DigiDepositSmearEnergy::smear);
0099       }
0100
0101       /// Create deposit mapping with updates on same cellIDs
0102       template <typename T> void
0103       smear(DigiContext& context, T& cont, work_t& /* work */, const predicate_t& predicate)  const  {
0104         auto& random = context.randomGenerator();
0105         std::size_t updated = 0UL;
0106
0107         for( auto& dep : cont )    {
0108           if ( predicate(dep) )   {
0109             CellID cell = dep.first;
0110             EnergyDeposit& depo = dep.second;
0111             double deposit = depo.deposit;
0112             double delta_E = 0e0;
0113             double energy  = deposit / dd4hep::GeV; // E in units of GeV
0114             double sigma_E_systematic   = m_systematic_resolution * energy;
0115             double sigma_E_intrin_fluct = m_intrinsic_fluctuation * std::sqrt(energy);
0116             double sigma_E_instrument   = m_instrumentation_resolution / dd4hep::GeV;
0117             double delta_ion = 0e0, num_pairs = 0e0;
0118             constexpr static double eps = std::numeric_limits<double>::epsilon();
0119             if ( sigma_E_systematic > eps )   {
0120               delta_E += sigma_E_systematic * random.gaussian(0e0, 1e0);
0121             }
0122             if ( sigma_E_intrin_fluct > eps )   {
0123               delta_E += sigma_E_intrin_fluct * random.gaussian(0e0, 1e0);
0124             }
0125             if ( sigma_E_instrument > eps )   {
0126               delta_E += sigma_E_instrument * random.gaussian(0e0, 1e0);
0127             }
0128             if ( m_ionization_fluctuation )   {
0129               num_pairs = energy / (m_pair_ionization_energy/dd4hep::GeV);
0130               delta_ion = energy * (random.poisson(num_pairs)/num_pairs);
0131               delta_E += delta_ion;
0132             }
0133             if ( dd4hep::isActivePrintLevel(outputLevel()) )   {
0134               print("%s+++ %016lX [GeV] E:%9.2e [%9.2e %9.2e] intrin_fluct:%9.2e systematic:%9.2e instrument:%9.2e ioni:%9.2e/%.0f",
0135                     context.event->id(), cell, energy, deposit/dd4hep::GeV, delta_E,
0136                     sigma_E_intrin_fluct, sigma_E_systematic, sigma_E_instrument, delta_ion, num_pairs);
0137             }
0138             /// delta_E is in GeV
0139             delta_E *= dd4hep::GeV;
0140             depo.depositError = delta_E;
0141             if ( m_monitor )  {
0142               m_monitor->energy_shift(dep, delta_E);
0143             }
0144             if ( m_modify_energy )  {
0145               depo.deposit = deposit + delta_E;
0146               depo.flag |= EnergyDeposit::ENERGY_SMEARED;
0147             }
0148             ++updated;
0149           }
0150         }
0151         info("%s+++ %-32s Smear energy: updated %6ld out of %6ld entries from mask: %04X",
0152              context.event->id(), cont.name.c_str(), updated, cont.size(), cont.key.mask());
0153       }
0154     };
0155
0156     /// Actor to only set energy error (as above, but with preset option
0157     /**
0158      *
0159      *  \author  M.Frank
0160      *  \version 1.0
0161      *  \ingroup DD4HEP_DIGITIZATION
0162      */
0163     class DigiDepositSetEnergyError : public DigiDepositSmearEnergy   {
0164     public:
0165       /// Standard constructor
0166       DigiDepositSetEnergyError(const DigiKernel& krnl, const std::string& nam)
0167         : DigiDepositSmearEnergy(krnl, nam)
0168       {
0169         m_modify_energy = false;
0170       }
0171     };
0172   }    // End namespace digi
0173 }      // End namespace dd4hep
0174
0175 #include <DDDigi/DigiFactories.h>
0176 DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDepositSmearEnergy)
0177 DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDepositSetEnergyError)