src/components/CalorimeterBirksCorr.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Wouter Deconinck, Whitney Armstrong, Sylvester Joosten, Jihee Kim
0003
0004 // Apply Birks Law to correct the energy deposit
0005 // It uses the contributions member in edm4hep::CalorimeterHit, so simulation must enable storeCalorimeterContributions
0006 //
0007 // Author: Chao Peng
0008 // Date: 09/29/2021
0009
0010 #include <algorithm>
0011 #include <cmath>
0012
0013 #include "Gaudi/Algorithm.h"
0014 #include "GaudiKernel/PhysicalConstants.h"
0015 #include "Gaudi/Property.h"
0016 #include "GaudiKernel/RndmGenerators.h"
0017
0018 #include <k4FWCore/DataHandle.h>
0019 #include "JugBase/IParticleSvc.h"
0020
0021 // Event Model related classes
0022 #include "edm4hep/SimCalorimeterHitCollection.h"
0023
0024
0025 using namespace Gaudi::Units;
0026
0027 namespace Jug::Digi {
0028
0029   /** Generic calorimeter hit digitiziation.
0030    *
0031    * \ingroup digi
0032    * \ingroup calorimetry
0033    */
0034   class CalorimeterBirksCorr : public Gaudi::Algorithm {
0035   public:
0036
0037     // digitization settings
0038     Gaudi::Property<double> m_birksConstant{this, "birksConstant", 0.126*mm/MeV};
0039
0040     mutable DataHandle<edm4hep::SimCalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
0041                                                                       this};
0042     mutable DataHandle<edm4hep::SimCalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0043                                                                        this};
0044
0045     SmartIF<IParticleSvc> m_pidSvc;
0046     // unitless conterpart of arguments
0047     double birksConstant{0};
0048
0049     //  ill-formed: using Gaudi::Algorithm::GaudiAlgorithm;
0050     CalorimeterBirksCorr(const std::string& name, ISvcLocator* svcLoc)
0051       : Gaudi::Algorithm(name, svcLoc)
0052     {
0053       declareProperty("inputHitCollection", m_inputHitCollection, "");
0054       declareProperty("outputHitCollection", m_outputHitCollection, "");
0055     }
0056
0057     StatusCode initialize() override
0058     {
0059       if (Gaudi::Algorithm::initialize().isFailure()) {
0060         return StatusCode::FAILURE;
0061       }
0062
0063       m_pidSvc = service("ParticleSvc");
0064       if (!m_pidSvc) {
0065         error() << "Unable to locate Particle Service. "
0066                 << "Make sure you have ParticleSvc in the configuration."
0067                 << endmsg;
0068         return StatusCode::FAILURE;
0069       }
0070
0071       // using juggler internal units (GeV, mm, radian, ns)
0072       birksConstant = m_birksConstant.value() / mm * GeV;
0073
0074       return StatusCode::SUCCESS;
0075     }
0076
0077     StatusCode execute(const EventContext&) const override
0078     {
0079       auto& ohits = *m_outputHitCollection.createAndPut();
0080       for (const auto& hit : *m_inputHitCollection.get()) {
0081         auto ohit = ohits->create(hit.getCellID(), hit.getEnergy(), hit.getPosition());
0082         double energy = 0.;
0083         for (const auto &c: hit.getContributions()) {
0084           ohit.addToContributions(c);
0085           const double charge = m_pidSvc->particle(c.getPDG()).charge;
0086           // some tolerance for precision
0087           if (std::abs(charge) > 1e-5) {
0088             // FIXME
0089             //energy += c.getEnergy() / (1. + c.getEnergy() / c.length * birksConstant);
0090             error() << "edm4hep::CaloHitContribution has no length field for Birks correction." << endmsg;
0091             return StatusCode::FAILURE;
0092           }
0093         }
0094         // replace energy deposit with Birks Law corrected value
0095         ohit.setEnergy(energy);
0096       }
0097       return StatusCode::SUCCESS;
0098     }
0099   };
0100   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0101   DECLARE_COMPONENT(CalorimeterBirksCorr)
0102
0103 } // namespace Jug::Digi