File indexing completed on 2023-03-17 11:24:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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 "G4SystemOfUnits.hh"
0027
0028
0029
0030
0031
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
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
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 }