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