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 "GaudiAlg/GaudiAlgorithm.h"
0014 #include "GaudiAlg/GaudiTool.h"
0015 #include "GaudiAlg/Transformer.h"
0016 #include "GaudiKernel/PhysicalConstants.h"
0017 #include "Gaudi/Property.h"
0018 #include "GaudiKernel/RndmGenerators.h"
0019
0020 #include <k4FWCore/DataHandle.h>
0021 #include "JugBase/IParticleSvc.h"
0022
0023 // Event Model related classes
0024 #include "edm4hep/SimCalorimeterHitCollection.h"
0025
0026
0027 using namespace Gaudi::Units;
0028
0029 namespace Jug::Digi {
0030
0031   /** Generic calorimeter hit digitiziation.
0032    *
0033    * \ingroup digi
0034    * \ingroup calorimetry
0035    */
0036   class CalorimeterBirksCorr : public GaudiAlgorithm {
0037   public:
0038
0039     // digitization settings
0040     Gaudi::Property<double> m_birksConstant{this, "birksConstant", 0.126*mm/MeV};
0041
0042     DataHandle<edm4hep::SimCalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
0043                                                                       this};
0044     DataHandle<edm4hep::SimCalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0045                                                                        this};
0046
0047     SmartIF<IParticleSvc> m_pidSvc;
0048     // unitless conterpart of arguments
0049     double birksConstant{0};
0050
0051     //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
0052     CalorimeterBirksCorr(const std::string& name, ISvcLocator* svcLoc)
0053       : GaudiAlgorithm(name, svcLoc)
0054     {
0055       declareProperty("inputHitCollection", m_inputHitCollection, "");
0056       declareProperty("outputHitCollection", m_outputHitCollection, "");
0057     }
0058
0059     StatusCode initialize() override
0060     {
0061       if (GaudiAlgorithm::initialize().isFailure()) {
0062         return StatusCode::FAILURE;
0063       }
0064
0065       m_pidSvc = service("ParticleSvc");
0066       if (!m_pidSvc) {
0067         error() << "Unable to locate Particle Service. "
0068                 << "Make sure you have ParticleSvc in the configuration."
0069                 << endmsg;
0070         return StatusCode::FAILURE;
0071       }
0072
0073       // using juggler internal units (GeV, mm, radian, ns)
0074       birksConstant = m_birksConstant.value() / mm * GeV;
0075
0076       return StatusCode::SUCCESS;
0077     }
0078
0079     StatusCode execute() override
0080     {
0081       auto& ohits = *m_outputHitCollection.createAndPut();
0082       for (const auto& hit : *m_inputHitCollection.get()) {
0083         auto ohit = ohits->create(hit.getCellID(), hit.getEnergy(), hit.getPosition());
0084         double energy = 0.;
0085         for (const auto &c: hit.getContributions()) {
0086           ohit.addToContributions(c);
0087           const double charge = m_pidSvc->particle(c.getPDG()).charge;
0088           // some tolerance for precision
0089           if (std::abs(charge) > 1e-5) {
0090             // FIXME
0091             //energy += c.getEnergy() / (1. + c.getEnergy() / c.length * birksConstant);
0092             error() << "edm4hep::CaloHitContribution has no length field for Birks correction." << endmsg;
0093             return StatusCode::FAILURE;
0094           }
0095         }
0096         // replace energy deposit with Birks Law corrected value
0097         ohit.setEnergy(energy);
0098       }
0099       return StatusCode::SUCCESS;
0100     }
0101   };
0102   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0103   DECLARE_COMPONENT(CalorimeterBirksCorr)
0104
0105 } // namespace Jug::Digi