Geometry/src/MaterialBudgetHcal.cc

0001 #include "Validation/Geometry/interface/MaterialBudgetHcal.h"
0002
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0008 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0009
0010 #include "G4Step.hh"
0011 #include "G4Track.hh"
0012
0013 #include <iostream>
0014 #include <memory>
0015
0016 MaterialBudgetHcal::MaterialBudgetHcal(const edm::ParameterSet& p) {
0017   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MaterialBudgetHcal");
0018   rMax_ = m_p.getUntrackedParameter<double>("RMax", 4.5) * CLHEP::m;
0019   zMax_ = m_p.getUntrackedParameter<double>("ZMax", 13.0) * CLHEP::m;
0020   fromdd4hep_ = m_p.getUntrackedParameter<bool>("Fromdd4hep", false);
0021   bool doHcal = m_p.getUntrackedParameter<bool>("DoHCAL", true);
0022   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcal initialized with rMax " << rMax_ << " mm and zMax " << zMax_
0023                                      << " mm doHcal is set to " << doHcal << " and Fromdd4hep to " << fromdd4hep_;
0024   if (doHcal) {
0025     theHistoHcal_ = std::make_unique<MaterialBudgetHcalHistos>(m_p);
0026     theHistoCastor_.reset(nullptr);
0027   } else {
0028     theHistoHcal_.reset(nullptr);
0029     theHistoCastor_ = std::make_unique<MaterialBudgetCastorHistos>(m_p);
0030   }
0031 }
0032
0033 void MaterialBudgetHcal::registerConsumes(edm::ConsumesCollector cc) {
0034   if (fromdd4hep_)
0035     cpvTokenDD4hep_ = cc.esConsumes<edm::Transition::BeginRun>();
0036   else
0037     cpvTokenDDD_ = cc.esConsumes<edm::Transition::BeginRun>();
0038   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcal: Initialize the token for CompactView";
0039 }
0040
0041 void MaterialBudgetHcal::beginRun(edm::EventSetup const& es) {
0042   //----- Check that selected volumes are indeed part of the geometry
0043   // Numbering From DDD
0044   if (fromdd4hep_) {
0045     const auto& pDD = es.getHandle(cpvTokenDD4hep_);
0046     if (theHistoHcal_)
0047       theHistoHcal_->fillBeginJob((*pDD));
0048   } else {
0049     const auto& pDD = es.getHandle(cpvTokenDDD_);
0050     if (theHistoHcal_)
0051       theHistoHcal_->fillBeginJob((*pDD));
0052   }
0053 }
0054
0055 void MaterialBudgetHcal::update(const BeginOfTrack* trk) {
0056   const G4Track* aTrack = (*trk)();  // recover G4 pointer if wanted
0057   if (theHistoHcal_)
0058     theHistoHcal_->fillStartTrack(aTrack);
0059   if (theHistoCastor_)
0060     theHistoCastor_->fillStartTrack(aTrack);
0061 }
0062
0063 void MaterialBudgetHcal::update(const G4Step* aStep) {
0064   //---------- each step
0065   if (theHistoHcal_)
0066     theHistoHcal_->fillPerStep(aStep);
0067   if (theHistoCastor_)
0068     theHistoCastor_->fillPerStep(aStep);
0069
0070   //----- Stop tracking after selected position
0071   if (stopAfter(aStep)) {
0072     G4Track* track = aStep->GetTrack();
0073     track->SetTrackStatus(fStopAndKill);
0074   }
0075 }
0076
0077 void MaterialBudgetHcal::update(const EndOfTrack* trk) {
0078   if (theHistoHcal_)
0079     theHistoHcal_->fillEndTrack();
0080   if (theHistoCastor_)
0081     theHistoCastor_->fillEndTrack();
0082 }
0083
0084 bool MaterialBudgetHcal::stopAfter(const G4Step* aStep) {
0085   G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0086   double rr = hitPoint.perp();
0087   double zz = std::abs(hitPoint.z());
0088
0089   if (rr > rMax_ || zz > zMax_) {
0090     edm::LogVerbatim("MaterialBudget") << " MaterialBudgetHcal::StopAfter R = " << rr << " and Z = " << zz;
0091     return true;
0092   } else {
0093     return false;
0094   }
0095 }