HGCalTestBeam/plugins/HGCalTBMBProducer.cc

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 //#define EDM_ML_DEBUG
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;                   // stop default
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)();  // recover G4 pointer if wanted
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);