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