File indexing completed on 2023-03-17 11:24:35
0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ServiceRegistry/interface/Service.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006
0007 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0008 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0009 #include "SimG4Core/Notification/interface/Observer.h"
0010 #include "SimG4Core/Watcher/interface/SimProducer.h"
0011 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingCalo.h"
0012
0013 #include "G4LogicalVolumeStore.hh"
0014 #include "G4Step.hh"
0015 #include "G4Track.hh"
0016
0017 #include <iostream>
0018 #include <string>
0019 #include <vector>
0020
0021
0022
0023 class HGCalTBMBProducer : public SimProducer,
0024 public Observer<const BeginOfTrack*>,
0025 public Observer<const G4Step*>,
0026 public Observer<const EndOfTrack*> {
0027 public:
0028 HGCalTBMBProducer(const edm::ParameterSet&);
0029 HGCalTBMBProducer(const HGCalTBMBProducer&) = delete;
0030 const HGCalTBMBProducer& operator=(const HGCalTBMBProducer&) = delete;
0031 ~HGCalTBMBProducer() override = default;
0032
0033 void produce(edm::Event&, const edm::EventSetup&) override;
0034
0035 private:
0036 void update(const BeginOfTrack*) override;
0037 void update(const G4Step*) override;
0038 void update(const EndOfTrack*) override;
0039
0040 bool stopAfter(const G4Step*);
0041 int findVolume(const G4VTouchable* touch, bool stop) const;
0042
0043 const edm::ParameterSet m_p;
0044 const std::vector<std::string> listNames_;
0045 const std::string stopName_;
0046 const double stopZ_;
0047 const unsigned int nList_;
0048 MaterialAccountingCaloCollection matcoll_;
0049 std::vector<double> radLen_, intLen_, stepLen_;
0050 };
0051
0052 HGCalTBMBProducer::HGCalTBMBProducer(const edm::ParameterSet& p)
0053 : m_p(p.getParameter<edm::ParameterSet>("HGCalTBMB")),
0054 listNames_(m_p.getParameter<std::vector<std::string> >("DetectorNames")),
0055 stopName_(m_p.getParameter<std::string>("StopName")),
0056 stopZ_(m_p.getParameter<double>("MaximumZ")),
0057 nList_(listNames_.size()) {
0058 edm::LogVerbatim("HGCSim") << "HGCalTBMBProducer initialized for " << nList_ << " volumes";
0059 for (unsigned int k = 0; k < nList_; ++k)
0060 edm::LogVerbatim("HGCSim") << " [" << k << "] " << listNames_[k];
0061 edm::LogVerbatim("HGCSim") << "Stop after " << stopZ_ << " or reaching volume " << stopName_;
0062
0063 produces<MaterialAccountingCaloCollection>("HGCalTBMB");
0064 }
0065
0066 void HGCalTBMBProducer::produce(edm::Event& e, const edm::EventSetup&) {
0067 std::unique_ptr<MaterialAccountingCaloCollection> hgc(new MaterialAccountingCaloCollection);
0068 for (auto const& mbc : matcoll_) {
0069 hgc->emplace_back(mbc);
0070 }
0071 e.put(std::move(hgc), "HGCalTBMB");
0072 }
0073
0074 void HGCalTBMBProducer::update(const BeginOfTrack* trk) {
0075 radLen_ = std::vector<double>(nList_ + 1, 0);
0076 intLen_ = std::vector<double>(nList_ + 1, 0);
0077 stepLen_ = std::vector<double>(nList_ + 1, 0);
0078
0079 #ifdef EDM_ML_DEBUG
0080 const G4Track* aTrack = (*trk)();
0081 const G4ThreeVector& mom = aTrack->GetMomentum();
0082 double theEnergy = aTrack->GetTotalEnergy();
0083 int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
0084 edm::LogVerbatim("HGCSim") << "HGCalTBMBProducer: Track " << aTrack->GetTrackID() << " Code " << theID << " Energy "
0085 << theEnergy / CLHEP::GeV << " GeV; Momentum " << mom;
0086 #endif
0087 }
0088
0089 void HGCalTBMBProducer::update(const G4Step* aStep) {
0090 G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
0091 double step = aStep->GetStepLength();
0092 double radl = material->GetRadlen();
0093 double intl = material->GetNuclearInterLength();
0094
0095 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0096 int indx = findVolume(touch, false);
0097
0098 if (indx >= 0) {
0099 stepLen_[indx] += step;
0100 radLen_[indx] += (step / radl);
0101 intLen_[indx] += (step / intl);
0102 }
0103 stepLen_[nList_] += step;
0104 radLen_[nList_] += (step / radl);
0105 intLen_[nList_] += (step / intl);
0106 #ifdef EDM_ML_DEBUG
0107 edm::LogVerbatim("HGCSim") << "HGCalTBMBProducer::Step in " << touch->GetVolume(0)->GetLogicalVolume()->GetName()
0108 << " Index " << indx << " Step " << step << " RadL " << step / radl << " IntL "
0109 << step / intl;
0110 #endif
0111
0112 if (stopAfter(aStep)) {
0113 G4Track* track = aStep->GetTrack();
0114 track->SetTrackStatus(fStopAndKill);
0115 }
0116 }
0117
0118 void HGCalTBMBProducer::update(const EndOfTrack* trk) {
0119 MaterialAccountingCalo matCalo;
0120 matCalo.m_eta = 0;
0121 matCalo.m_phi = 0;
0122 for (unsigned int ii = 0; ii <= nList_; ++ii) {
0123 matCalo.m_stepLen.emplace_back(stepLen_[ii]);
0124 matCalo.m_radLen.emplace_back(radLen_[ii]);
0125 matCalo.m_intLen.emplace_back(intLen_[ii]);
0126 #ifdef EDM_ML_DEBUG
0127 std::string name("Total");
0128 if (ii < nList_)
0129 name = listNames_[ii];
0130 edm::LogVerbatim("HGCSim") << "HGCalTBMBProducer::Volume[" << ii << "]: " << name << " == Step " << stepLen_[ii]
0131 << " RadL " << radLen_[ii] << " IntL " << intLen_[ii];
0132 #endif
0133 }
0134 matcoll_.emplace_back(matCalo);
0135 }
0136
0137 bool HGCalTBMBProducer::stopAfter(const G4Step* aStep) {
0138 bool flag(false);
0139 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0140 G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0141 if (aStep->GetPostStepPoint() != nullptr)
0142 hitPoint = aStep->GetPostStepPoint()->GetPosition();
0143 double zz = hitPoint.z();
0144
0145 if ((findVolume(touch, true) == 0) || (zz > stopZ_))
0146 flag = true;
0147 #ifdef EDM_ML_DEBUG
0148 edm::LogVerbatim("HGCSim") << " HGCalTBMBProducer::Name " << touch->GetVolume(0)->GetName() << " z " << zz << " Flag"
0149 << flag;
0150 #endif
0151 return flag;
0152 }
0153
0154 int HGCalTBMBProducer::findVolume(const G4VTouchable* touch, bool stop) const {
0155 int ivol = -1;
0156 int level = (touch->GetHistoryDepth()) + 1;
0157 for (int ii = 0; ii < level; ii++) {
0158 std::string name = touch->GetVolume(ii)->GetName();
0159 if (stop) {
0160 if (strcmp(name.c_str(), stopName_.c_str()) == 0)
0161 ivol = 0;
0162 } else {
0163 for (unsigned int k = 0; k < nList_; ++k) {
0164 if (strcmp(name.c_str(), listNames_[k].c_str()) == 0) {
0165 ivol = k;
0166 break;
0167 }
0168 }
0169 }
0170 if (ivol >= 0)
0171 break;
0172 }
0173 return ivol;
0174 }
0175
0176 #include "FWCore/Framework/interface/MakerMacros.h"
0177 #include "FWCore/PluginManager/interface/ModuleDef.h"
0178 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0179
0180 DEFINE_SIMWATCHER(HGCalTBMBProducer);