Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-24 02:22:26

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/EventSetup.h"
0003 #include "FWCore/Framework/interface/ESTransientHandle.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 
0009 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0010 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0011 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0012 #include "SimG4Core/Watcher/interface/SimProducer.h"
0013 #include "SimG4Core/Notification/interface/Observer.h"
0014 
0015 #include "DataFormats/Math/interface/GeantUnits.h"
0016 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0017 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingCalo.h"
0018 
0019 #include "DetectorDescription/Core/interface/DDCompactView.h"
0020 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0021 #include "DetectorDescription/Core/interface/DDFilter.h"
0022 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0023 #include "DetectorDescription/Core/interface/DDSplit.h"
0024 #include "DetectorDescription/Core/interface/DDValue.h"
0025 #include "DetectorDescription/Core/interface/DDsvalues.h"
0026 #include "DetectorDescription/DDCMS/interface/DDCompactView.h"
0027 #include "DetectorDescription/DDCMS/interface/DDFilteredView.h"
0028 
0029 #include <CLHEP/Vector/LorentzVector.h>
0030 
0031 #include "DD4hep/Filter.h"
0032 
0033 #include "G4Step.hh"
0034 #include "G4Track.hh"
0035 
0036 #include <iostream>
0037 #include <memory>
0038 #include <string>
0039 #include <vector>
0040 
0041 using namespace geant_units::operators;
0042 
0043 class MaterialBudgetHcalProducer : public SimProducer,
0044                                    public Observer<const BeginOfEvent*>,
0045                                    public Observer<const BeginOfTrack*>,
0046                                    public Observer<const G4Step*>,
0047                                    public Observer<const EndOfTrack*> {
0048 public:
0049   MaterialBudgetHcalProducer(const edm::ParameterSet&);
0050   ~MaterialBudgetHcalProducer() override = default;
0051 
0052   MaterialBudgetHcalProducer(const MaterialBudgetHcalProducer&) = delete;                   // stop default
0053   const MaterialBudgetHcalProducer& operator=(const MaterialBudgetHcalProducer&) = delete;  // stop default
0054 
0055   void registerConsumes(edm::ConsumesCollector) override;
0056   void produce(edm::Event&, const edm::EventSetup&) override;
0057   void beginRun(edm::EventSetup const&) override;
0058 
0059 private:
0060   void update(const BeginOfEvent*) override;
0061   void update(const BeginOfTrack*) override;
0062   void update(const G4Step*) override;
0063   void update(const EndOfTrack*) override;
0064 
0065   bool stopAfter(const G4Step*);
0066   std::vector<std::string> getNames(DDFilteredView& fv);
0067   std::vector<std::string> getNames(cms::DDFilteredView& fv);
0068   std::vector<double> getDDDArray(const std::string& str, const DDsvalues_type& sv);
0069   bool isSensitive(const std::string&);
0070   bool isItHF(const G4VTouchable*);
0071   bool isItEC(const std::string&);
0072 
0073   static const int maxSet_ = 25, maxSet2_ = 9;
0074   double rMax_, zMax_, etaMinP_, etaMaxP_;
0075   bool printSum_, fromdd4hep_;
0076   std::string name_;
0077 
0078   edm::ESGetToken<DDCompactView, IdealGeometryRecord> cpvTokenDDD_;
0079   edm::ESGetToken<cms::DDCompactView, IdealGeometryRecord> cpvTokenDD4hep_;
0080 
0081   std::vector<std::string> sensitives_, hfNames_, sensitiveEC_;
0082   std::vector<int> hfLevels_;
0083   MaterialAccountingCaloCollection matcoll_;
0084   double stepLens_[maxSet_], radLens_[maxSet_], intLens_[maxSet_];
0085   std::vector<std::string> matList_;
0086   std::vector<double> stepLength_, radLength_, intLength_;
0087   int id_, layer_, steps_;
0088   double radLen_, intLen_, stepLen_;
0089   double eta_, phi_;
0090   int nlayHB_, nlayHE_, nlayHO_, nlayHF_;
0091 };
0092 
0093 MaterialBudgetHcalProducer::MaterialBudgetHcalProducer(const edm::ParameterSet& p) {
0094   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MaterialBudgetHcalProducer");
0095   rMax_ = m_p.getUntrackedParameter<double>("RMax", 4.5) * CLHEP::m;
0096   zMax_ = m_p.getUntrackedParameter<double>("ZMax", 13.0) * CLHEP::m;
0097   etaMinP_ = m_p.getUntrackedParameter<double>("EtaMinP", 5.2);
0098   etaMaxP_ = m_p.getUntrackedParameter<double>("EtaMaxP", 0.0);
0099   fromdd4hep_ = m_p.getUntrackedParameter<bool>("Fromdd4hep", false);
0100   printSum_ = m_p.getUntrackedParameter<bool>("PrintSummary", false);
0101   name_ = m_p.getUntrackedParameter<std::string>("Name", "Hcal");
0102   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalProducer initialized with rMax " << rMax_ << " mm and zMax "
0103                                      << zMax_ << " mm printSummary is set to " << printSum_ << " and Fromdd4hep to "
0104                                      << fromdd4hep_;
0105 
0106   produces<MaterialAccountingCaloCollection>(Form("%sMatBCalo", name_.c_str()));
0107 }
0108 
0109 void MaterialBudgetHcalProducer::registerConsumes(edm::ConsumesCollector cc) {
0110   if (fromdd4hep_)
0111     cpvTokenDD4hep_ = cc.esConsumes<edm::Transition::BeginRun>();
0112   else
0113     cpvTokenDDD_ = cc.esConsumes<edm::Transition::BeginRun>();
0114   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalProducer: Initialize the token for CompactView";
0115 }
0116 
0117 void MaterialBudgetHcalProducer::produce(edm::Event& e, const edm::EventSetup&) {
0118   std::unique_ptr<MaterialAccountingCaloCollection> hgc(new MaterialAccountingCaloCollection);
0119   for (auto const& mbc : matcoll_) {
0120     hgc->emplace_back(mbc);
0121   }
0122   e.put(std::move(hgc), Form("%sMatBCalo", name_.c_str()));
0123 }
0124 
0125 void MaterialBudgetHcalProducer::beginRun(edm::EventSetup const& es) {
0126   //----- Check that selected volumes are indeed part of the geometry
0127   // Numbering From DDD
0128   if (fromdd4hep_) {
0129     const cms::DDCompactView cpv = es.getData(cpvTokenDD4hep_);
0130     constexpr int32_t addLevel = 1;
0131     std::string attribute = "ReadOutName";
0132     std::string value = "HcalHits";
0133     const cms::DDFilter filter1(attribute, value);
0134     cms::DDFilteredView fv1(cpv, filter1);
0135     std::vector<std::string> names = getNames(fv1);
0136     for (auto& name : names) {
0137       std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
0138       if (std::find(sensitives_.begin(), sensitives_.end(), namx) == sensitives_.end())
0139         sensitives_.emplace_back(namx);
0140     }
0141     edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalProducer: Names to be tested for " << attribute
0142                                            << " = " << value << " has " << sensitives_.size() << " elements";
0143     for (unsigned int i = 0; i < sensitives_.size(); i++)
0144       edm::LogVerbatim("MaterialBudgetFull")
0145           << "MaterialBudgetHcalProducer: sensitives[" << i << "] = " << sensitives_[i];
0146     attribute = "Volume";
0147     value = "HF";
0148     const cms::DDFilter filter2(attribute, value);
0149     cms::DDFilteredView fv2(cpv, filter2);
0150     std::vector<int> temp = fv2.get<std::vector<int> >("hf", "Levels");
0151     hfNames_ = getNames(fv2);
0152     edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalProducer: Names to be tested for " << attribute
0153                                            << " = " << value << " has " << hfNames_.size() << " elements";
0154     for (unsigned int i = 0; i < hfNames_.size(); i++) {
0155       hfLevels_.push_back(temp[i] + addLevel);
0156       edm::LogVerbatim("MaterialBudgetFull")
0157           << "MaterialBudgetHcalProducer:  HF[" << i << "] = " << hfNames_[i] << " at level " << hfLevels_[i];
0158     }
0159 
0160     const std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
0161     attribute = "ReadOutName";
0162     for (int k = 0; k < 2; k++) {
0163       value = ecalRO[k];
0164       const cms::DDFilter filter(attribute, value);
0165       cms::DDFilteredView fv(cpv, filter);
0166       std::vector<std::string> senstmp = getNames(fv);
0167       edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalProducer: Names to be tested for " << attribute
0168                                              << " = " << value << " has " << senstmp.size() << " elements";
0169       for (unsigned int i = 0; i < senstmp.size(); i++) {
0170         std::string name = senstmp[i].substr(0, 4);
0171         if (std::find(sensitiveEC_.begin(), sensitiveEC_.end(), name) == sensitiveEC_.end())
0172           sensitiveEC_.push_back(name);
0173       }
0174     }
0175     for (unsigned int i = 0; i < sensitiveEC_.size(); i++)
0176       edm::LogVerbatim("MaterialBudgetFull")
0177           << "MaterialBudgetHcalProducer:sensitiveEC[" << i << "] = " << sensitiveEC_[i];
0178   } else {  // if not from dd4hep --> ddd
0179     const DDCompactView& cpv = es.getData(cpvTokenDDD_);
0180     constexpr int32_t addLevel = 0;
0181     std::string attribute = "ReadOutName";
0182     std::string value = "HcalHits";
0183     DDSpecificsMatchesValueFilter filter1{DDValue(attribute, value, 0)};
0184     DDFilteredView fv1(cpv, filter1);
0185     std::vector<std::string> names = getNames(fv1);
0186     for (auto& name : names) {
0187       std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
0188       if (std::find(sensitives_.begin(), sensitives_.end(), namx) == sensitives_.end())
0189         sensitives_.emplace_back(namx);
0190     }
0191     edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalProducer: Names to be tested for " << attribute
0192                                            << " = " << value << " has " << sensitives_.size() << " elements";
0193     for (unsigned int i = 0; i < sensitives_.size(); i++)
0194       edm::LogVerbatim("MaterialBudgetFull")
0195           << "MaterialBudgetHcalProducer: sensitives[" << i << "] = " << sensitives_[i];
0196     attribute = "Volume";
0197     value = "HF";
0198     DDSpecificsMatchesValueFilter filter2{DDValue(attribute, value, 0)};
0199     DDFilteredView fv2(cpv, filter2);
0200     hfNames_ = getNames(fv2);
0201     fv2.firstChild();
0202     DDsvalues_type sv(fv2.mergedSpecifics());
0203     std::vector<double> temp = getDDDArray("Levels", sv);
0204     edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalProducer: Names to be tested for " << attribute
0205                                            << " = " << value << " has " << hfNames_.size() << " elements";
0206     for (unsigned int i = 0; i < hfNames_.size(); i++) {
0207       int level = static_cast<int>(temp[i]);
0208       hfLevels_.push_back(level + addLevel);
0209       edm::LogVerbatim("MaterialBudgetFull")
0210           << "MaterialBudgetHcalProducer:  HF[" << i << "] = " << hfNames_[i] << " at level " << hfLevels_[i];
0211     }
0212 
0213     const std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
0214     attribute = "ReadOutName";
0215     for (int k = 0; k < 2; k++) {
0216       value = ecalRO[k];
0217       DDSpecificsMatchesValueFilter filter3{DDValue(attribute, value, 0)};
0218       DDFilteredView fv3(cpv, filter3);
0219       std::vector<std::string> senstmp = getNames(fv3);
0220       edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalProducer: Names to be tested for " << attribute
0221                                              << " = " << value << " has " << senstmp.size() << " elements";
0222       for (unsigned int i = 0; i < senstmp.size(); i++) {
0223         std::string name = senstmp[i].substr(0, 4);
0224         if (std::find(sensitiveEC_.begin(), sensitiveEC_.end(), name) == sensitiveEC_.end())
0225           sensitiveEC_.push_back(name);
0226       }
0227     }
0228     for (unsigned int i = 0; i < sensitiveEC_.size(); i++)
0229       edm::LogVerbatim("MaterialBudgetFull")
0230           << "MaterialBudgetHcalProducer:sensitiveEC[" << i << "] = " << sensitiveEC_[i];
0231   }
0232 }
0233 
0234 void MaterialBudgetHcalProducer::update(const BeginOfEvent*) { matcoll_.clear(); }
0235 
0236 void MaterialBudgetHcalProducer::update(const BeginOfTrack* trk) {
0237   const G4Track* aTrack = (*trk)();  // recover G4 pointer if wanted
0238 
0239   id_ = layer_ = steps_ = 0;
0240   radLen_ = intLen_ = stepLen_ = 0;
0241   nlayHB_ = nlayHE_ = nlayHF_ = nlayHO_ = 0;
0242   for (int i = 0; i < maxSet_; ++i)
0243     stepLens_[i] = radLens_[i] = intLens_[i] = 0;
0244 
0245   const G4ThreeVector& dir = aTrack->GetMomentum();
0246   if (dir.theta() != 0) {
0247     eta_ = dir.eta();
0248   } else {
0249     eta_ = -99;
0250   }
0251   phi_ = dir.phi();
0252   edm::LogVerbatim("MaterialBudget") << "Start track with (eta, phi) = (" << eta_ << ", " << phi_ << " Energy "
0253                                      << aTrack->GetTotalEnergy() << " and ID "
0254                                      << (aTrack->GetDefinition()->GetPDGEncoding());
0255 
0256   if (printSum_) {
0257     matList_.clear();
0258     stepLength_.clear();
0259     radLength_.clear();
0260     intLength_.clear();
0261   }
0262 }
0263 
0264 void MaterialBudgetHcalProducer::update(const G4Step* aStep) {
0265   //---------- each step
0266   G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
0267   double step = aStep->GetStepLength();
0268   double radl = material->GetRadlen();
0269   double intl = material->GetNuclearInterLength();
0270   double density = convertUnitsTo(1._g_per_cm3, material->GetDensity());
0271 
0272   int idOld = id_;
0273   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0274   std::string name = (static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(0)->GetName())));
0275   std::string matName = (static_cast<std::string>(dd4hep::dd::noNamespace(material->GetName())));
0276   if (printSum_) {
0277     bool found = false;
0278     for (unsigned int ii = 0; ii < matList_.size(); ii++) {
0279       if (matList_[ii] == matName) {
0280         stepLength_[ii] += step;
0281         radLength_[ii] += (step / radl);
0282         intLength_[ii] += (step / intl);
0283         found = true;
0284         break;
0285       }
0286     }
0287     if (!found) {
0288       matList_.push_back(matName);
0289       stepLength_.push_back(step);
0290       radLength_.push_back(step / radl);
0291       intLength_.push_back(step / intl);
0292     }
0293     if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
0294       edm::LogVerbatim("MaterialBudget") << "Volume " << name << " id " << id_ << ":" << idOld << " Step " << step
0295                                          << " Material " << matName << " Old Length " << stepLen_ << " X0 "
0296                                          << step / radl << ":" << radLen_ << " Lambda " << step / intl << ":"
0297                                          << intLen_;
0298   } else {
0299     if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
0300       edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalProducer: Step at " << name << " id " << id_ << ":"
0301                                          << idOld << " Length " << step << " in " << matName << " of density "
0302                                          << density << " g/cc; Radiation Length " << radl << " mm; Interaction Length "
0303                                          << intl << " mm\n                          Position "
0304                                          << aStep->GetPreStepPoint()->GetPosition() << " Cylindrical R "
0305                                          << aStep->GetPreStepPoint()->GetPosition().perp() << " Length (so far) "
0306                                          << stepLen_ << " L/X0 " << step / radl << "/" << radLen_ << " L/Lambda "
0307                                          << step / intl << "/" << intLen_;
0308   }
0309 
0310   int det = 0, lay = 0;
0311   double abseta = std::abs(eta_);
0312   edm::LogVerbatim("MaterialBudgetFull") << "Volume " << name << ":" << matName << " EC:Sensitive:HF " << isItEC(name)
0313                                          << ":" << isSensitive(name) << ":" << isItHF(touch) << " Eta " << abseta
0314                                          << " HC " << ((touch->GetReplicaNumber(1)) / 1000) << ":"
0315                                          << ((touch->GetReplicaNumber(0) / 10) % 100 + 3) << " X0 "
0316                                          << (radLen_ + (step / radl)) << " Lambda " << (intLen_ + (step / intl));
0317 
0318   if (isItEC(name)) {
0319     det = 1;
0320     lay = 1;
0321   } else {
0322     if (isSensitive(name)) {
0323       if (isItHF(touch)) {
0324         det = 5;
0325         lay = 21;
0326         if (lay != layer_)
0327           ++nlayHF_;
0328       } else {
0329         det = (touch->GetReplicaNumber(1)) / 1000;
0330         lay = (touch->GetReplicaNumber(0) / 10) % 100 + 3;
0331         if (det == 4) {
0332           if (abseta < 1.479)
0333             lay = layer_ + 1;
0334           else
0335             lay--;
0336           if (lay < 3)
0337             lay = 3;
0338           if (lay == layer_)
0339             lay++;
0340           if (lay > 20)
0341             lay = 20;
0342           if (lay != layer_)
0343             ++nlayHE_;
0344         } else if (lay != layer_) {
0345           if (lay >= 20)
0346             ++nlayHO_;
0347           else
0348             ++nlayHB_;
0349         }
0350       }
0351       edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalProducer: Det " << det << " Layer " << lay << " Eta "
0352                                              << eta_ << " Phi " << convertRadToDeg(phi_);
0353     } else if (layer_ == 1) {
0354       det = -1;
0355       lay = 2;
0356     }
0357   }
0358   if (det != 0) {
0359     if (lay != layer_) {
0360       id_ = lay;
0361       layer_ = lay;
0362     }
0363   }
0364 
0365   if (id_ > idOld) {
0366     if ((abseta >= etaMinP_) && (abseta <= etaMaxP_))
0367       edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalProducer: Step at " << name << " calls filHisto with "
0368                                          << (id_ - 1);
0369     stepLens_[id_ - 1] = stepLen_;
0370     radLens_[id_ - 1] = radLen_;
0371     intLens_[id_ - 1] = intLen_;
0372   }
0373 
0374   stepLen_ += step;
0375   radLen_ += (step / radl);
0376   intLen_ += (step / intl);
0377   if (id_ == 21) {
0378     if (!isItHF(aStep->GetPostStepPoint()->GetTouchable())) {
0379       if ((abseta >= etaMinP_) && (abseta <= etaMaxP_))
0380         edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalProducer: After HF in " << name << ":"
0381                                            << static_cast<std::string>(dd4hep::dd::noNamespace(
0382                                                   aStep->GetPostStepPoint()->GetTouchable()->GetVolume(0)->GetName()))
0383                                            << " calls fillHisto with " << id_;
0384       stepLens_[idOld] = stepLen_;
0385       radLens_[idOld] = radLen_;
0386       intLens_[idOld] = intLen_;
0387       ++id_;
0388       layer_ = 0;
0389     }
0390   }
0391 
0392   //----- Stop tracking after selected position
0393   if (stopAfter(aStep)) {
0394     G4Track* track = aStep->GetTrack();
0395     track->SetTrackStatus(fStopAndKill);
0396   }
0397 }
0398 
0399 void MaterialBudgetHcalProducer::update(const EndOfTrack* trk) {
0400   MaterialAccountingCalo matCalo;
0401   matCalo.m_eta = eta_;
0402   matCalo.m_phi = phi_;
0403   for (int k = 0; k < maxSet_; ++k) {
0404     matCalo.m_stepLen.emplace_back(stepLens_[k]);
0405     matCalo.m_radLen.emplace_back(radLens_[k]);
0406     matCalo.m_intLen.emplace_back(intLens_[k]);
0407   }
0408   matCalo.m_layers.emplace_back(nlayHB_);
0409   matCalo.m_layers.emplace_back(nlayHE_);
0410   matCalo.m_layers.emplace_back(nlayHO_);
0411   matCalo.m_layers.emplace_back(nlayHF_);
0412   matcoll_.emplace_back(matCalo);
0413 }
0414 
0415 bool MaterialBudgetHcalProducer::stopAfter(const G4Step* aStep) {
0416   G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0417   double rr = hitPoint.perp();
0418   double zz = std::abs(hitPoint.z());
0419 
0420   if (rr > rMax_ || zz > zMax_) {
0421     edm::LogVerbatim("MaterialBudget") << " MaterialBudgetHcalProducer::StopAfter R = " << rr << " and Z = " << zz;
0422     return true;
0423   } else {
0424     return false;
0425   }
0426 }
0427 
0428 std::vector<std::string> MaterialBudgetHcalProducer::getNames(DDFilteredView& fv) {
0429   std::vector<std::string> tmp;
0430   bool dodet = fv.firstChild();
0431   while (dodet) {
0432     const DDLogicalPart& log = fv.logicalPart();
0433     std::string namx = log.name().name();
0434     if (std::find(tmp.begin(), tmp.end(), namx) == tmp.end())
0435       tmp.push_back(namx);
0436     dodet = fv.next();
0437   }
0438   return tmp;
0439 }
0440 
0441 std::vector<std::string> MaterialBudgetHcalProducer::getNames(cms::DDFilteredView& fv) {
0442   std::vector<std::string> tmp;
0443   const std::vector<std::string> notIn = {
0444       "CALO", "HCal", "MBBTL", "MBBTR", "MBBTC", "MBAT", "MBBT_R1M", "MBBT_R1P", "MBBT_R1MX", "MBBT_R1PX", "VCAL"};
0445   while (fv.firstChild()) {
0446     const std::string n{fv.name().data(), fv.name().size()};
0447     if (std::find(notIn.begin(), notIn.end(), n) == notIn.end()) {
0448       std::string::size_type pos = n.find(':');
0449       const std::string namx = (pos == std::string::npos) ? n : std::string(n, pos + 1, n.size() - 1);
0450       if (std::find(tmp.begin(), tmp.end(), namx) == tmp.end())
0451         tmp.push_back(namx);
0452     }
0453   }
0454   return tmp;
0455 }
0456 
0457 std::vector<double> MaterialBudgetHcalProducer::getDDDArray(const std::string& str, const DDsvalues_type& sv) {
0458   edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalProducer:getDDDArray called for " << str;
0459   DDValue value(str);
0460   if (DDfetch(&sv, value)) {
0461     edm::LogVerbatim("MaterialBudgetFull") << value;
0462     const std::vector<double>& fvec = value.doubles();
0463     int nval = fvec.size();
0464     if (nval < 1) {
0465       throw cms::Exception("MaterialBudgetHcalProducer") << "nval = " << nval << " < 1 for array " << str << "\n";
0466     }
0467 
0468     return fvec;
0469   } else {
0470     throw cms::Exception("MaterialBudgetHcalProducer") << "cannot get array " << str << "\n";
0471   }
0472 }
0473 
0474 bool MaterialBudgetHcalProducer::isSensitive(const std::string& name) {
0475   std::vector<std::string>::const_iterator it = sensitives_.begin();
0476   std::vector<std::string>::const_iterator itEnd = sensitives_.end();
0477   std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
0478   for (; it != itEnd; ++it)
0479     if (namx == *it)
0480       return true;
0481   return false;
0482 }
0483 
0484 bool MaterialBudgetHcalProducer::isItHF(const G4VTouchable* touch) {
0485   int levels = ((touch->GetHistoryDepth()) + 1);
0486   for (unsigned int it = 0; it < hfNames_.size(); it++) {
0487     if (levels >= hfLevels_[it]) {
0488       std::string name =
0489           (static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(levels - hfLevels_[it])->GetName())))
0490               .substr(0, 4);
0491       if (name == hfNames_[it]) {
0492         return true;
0493       }
0494     }
0495   }
0496   return false;
0497 }
0498 
0499 bool MaterialBudgetHcalProducer::isItEC(const std::string& name) {
0500   std::vector<std::string>::const_iterator it = sensitiveEC_.begin();
0501   std::vector<std::string>::const_iterator itEnd = sensitiveEC_.end();
0502   for (; it != itEnd; ++it)
0503     if (name.substr(0, 4) == *it)
0504       return true;
0505   return false;
0506 }
0507 
0508 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0509 #include "FWCore/PluginManager/interface/ModuleDef.h"
0510 
0511 DEFINE_SIMWATCHER(MaterialBudgetHcalProducer);