HcalTestBeam/src/HcalTB02SD.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     HcalTestBeam
0004 // Class  :     HcalTB02SD
0005 //
0006 // Implementation:
0007 //     Sensitive Detector class for Hcal Test Beam 2002 detectors
0008 //
0009 // Original Author:
0010 //         Created:  Sun 21 10:14:34 CEST 2006
0011 //
0012
0013 // system include files
0014
0015 // user include files
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02SD.h"
0019 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02HcalNumberingScheme.h"
0020 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02XtalNumberingScheme.h"
0021
0022 #include "G4Step.hh"
0023 #include "G4Track.hh"
0024 #include "G4VProcess.hh"
0025
0026 #include <CLHEP/Units/SystemOfUnits.h>
0027
0028 //#define EDM_ML_DEBUG
0029
0030 //
0031 // constructors and destructor
0032 //
0033
0034 HcalTB02SD::HcalTB02SD(const std::string& name,
0035                        const HcalTB02Parameters* es,
0036                        const SensitiveDetectorCatalog& clg,
0037                        edm::ParameterSet const& p,
0038                        const SimTrackManager* manager)
0039     : CaloSD(name, clg, p, manager), hcalTB02Parameters_(es) {
0040   numberingScheme_.reset(nullptr);
0041   edm::ParameterSet m_SD = p.getParameter<edm::ParameterSet>("HcalTB02SD");
0042   useBirk_ = m_SD.getUntrackedParameter<bool>("UseBirkLaw", false);
0043   birk1_ = m_SD.getUntrackedParameter<double>("BirkC1", 0.013) * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
0044   birk2_ = m_SD.getUntrackedParameter<double>("BirkC2", 0.0568);
0045   birk3_ = m_SD.getUntrackedParameter<double>("BirkC3", 1.75);
0046   useWeight_ = true;
0047
0048   HcalTB02NumberingScheme* scheme = nullptr;
0049   if (name == "EcalHitsEB") {
0050     scheme = dynamic_cast<HcalTB02NumberingScheme*>(new HcalTB02XtalNumberingScheme());
0051     useBirk_ = false;
0052   } else if (name == "HcalHits") {
0053     scheme = dynamic_cast<HcalTB02NumberingScheme*>(new HcalTB02HcalNumberingScheme());
0054     useWeight_ = false;
0055   } else {
0056     edm::LogWarning("HcalTBSim") << "HcalTB02SD: ReadoutName " << name << " not supported\n";
0057   }
0058
0059   if (scheme)
0060     setNumberingScheme(scheme);
0061 #ifdef EDM_ML_DEBUG
0062   edm::LogVerbatim("HcalTBSim") << "***************************************************\n"
0063                                 << "*                                                 *\n"
0064                                 << "* Constructing a HcalTB02SD  with name " << GetName() << "\n"
0065                                 << "*                                                 *\n"
0066                                 << "***************************************************";
0067   edm::LogVerbatim("HcalTBSim") << "HcalTB02SD:: Use of Birks law is set to      " << useBirk_
0068                                 << "        with three constants kB = " << birk1_ << ", C1 = " << birk2_
0069                                 << ", C2 = " << birk3_;
0070 #endif
0071 }
0072
0073 HcalTB02SD::~HcalTB02SD() {}
0074
0075 //
0076 // member functions
0077 //
0078
0079 double HcalTB02SD::getEnergyDeposit(const G4Step* aStep) {
0080   auto const preStepPoint = aStep->GetPreStepPoint();
0081   std::string nameVolume = static_cast<std::string>(preStepPoint->GetPhysicalVolume()->GetName());
0082
0083   // take into account light collection curve for crystals
0084   double weight = 1.;
0085   if (useWeight_)
0086     weight *= curve_LY(nameVolume, preStepPoint);
0087   if (useBirk_)
0088     weight *= getAttenuation(aStep, birk1_, birk2_, birk3_);
0089   double edep = aStep->GetTotalEnergyDeposit() * weight;
0090 #ifdef EDM_ML_DEBUG
0091   edm::LogVerbatim("HcalTBSim") << "HcalTB02SD:: " << nameVolume << " Light Collection Efficiency " << weight
0092                                 << " Weighted Energy Deposit " << edep / CLHEP::MeV << " MeV";
0093 #endif
0094   return edep;
0095 }
0096
0097 uint32_t HcalTB02SD::setDetUnitId(const G4Step* aStep) {
0098   return (numberingScheme_ == nullptr ? 0 : (uint32_t)(numberingScheme_->getUnitID(aStep)));
0099 }
0100
0101 void HcalTB02SD::setNumberingScheme(HcalTB02NumberingScheme* scheme) {
0102   if (scheme != nullptr) {
0103     edm::LogVerbatim("HcalTBSim") << "HcalTB02SD: updates numbering scheme for " << GetName();
0104     numberingScheme_.reset(scheme);
0105   }
0106 }
0107
0108 double HcalTB02SD::curve_LY(const std::string& nameVolume, const G4StepPoint* stepPoint) {
0109   double weight = 1.;
0110   G4ThreeVector localPoint = setToLocal(stepPoint->GetPosition(), stepPoint->GetTouchable());
0111   double crlength = crystalLength(nameVolume);
0112   double dapd = 0.5 * crlength - localPoint.z();
0113   if (dapd >= -0.1 || dapd <= crlength + 0.1) {
0114     if (dapd <= 100.)
0115       weight = 1.05 - dapd * 0.0005;
0116   } else {
0117     edm::LogWarning("HcalTBSim") << "HcalTB02SD: light coll curve : wrong "
0118                                  << "distance to APD " << dapd << " crlength = " << crlength
0119                                  << " crystal name = " << nameVolume << " z of localPoint = " << localPoint.z()
0120                                  << " take weight = " << weight;
0121   }
0122 #ifdef EDM_ML_DEBUG
0123   edm::LogVerbatim("HcalTBSim") << "HcalTB02SD, light coll curve : " << dapd << " crlength = " << crlength
0124                                 << " crystal name = " << nameVolume << " z of localPoint = " << localPoint.z()
0125                                 << " take weight = " << weight;
0126 #endif
0127   return weight;
0128 }
0129
0130 double HcalTB02SD::crystalLength(const std::string& name) {
0131   double length = 230.;
0132   std::map<std::string, double>::const_iterator it = hcalTB02Parameters_->lengthMap_.find(name);
0133   if (it != hcalTB02Parameters_->lengthMap_.end())
0134     length = it->second;
0135   return length;
0136 }