File indexing completed on 2021-09-13 05:55:08
0001 #include <iostream> // FIXME: switch to MessagLogger & friends
0002 #include <vector>
0003 #include <string>
0004 #include <cassert>
0005 #include <exception>
0006 #include <tuple>
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0015 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0016 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0017 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0018
0019 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0020
0021
0022 #include "G4Step.hh"
0023 #include "G4Track.hh"
0024 #include "G4VSolid.hh"
0025 #include "G4LogicalVolumeStore.hh"
0026 #include "G4TouchableHistory.hh"
0027 #include "G4VPhysicalVolume.hh"
0028 #include "G4AffineTransform.hh"
0029
0030 #include "TrackingMaterialProducer.h"
0031
0032
0033
0034
0035 #define DEBUG_G4_VOLUMES
0036
0037 using namespace CLHEP;
0038 using edm::LogInfo;
0039
0040
0041 static const G4LogicalVolume* GetVolume(const std::string& name) {
0042 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0043
0044 #ifdef DEBUG_G4_VOLUMES
0045 for (G4LogicalVolumeStore::const_iterator volume = lvs->begin(); volume != lvs->end(); ++volume)
0046 LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: G4 registered volumes " << (*volume)->GetName()
0047 << std::endl;
0048 #endif
0049
0050 for (G4LogicalVolumeStore::const_iterator volume = lvs->begin(); volume != lvs->end(); ++volume) {
0051 if ((const std::string&)(*volume)->GetName() == name)
0052 return (*volume);
0053 }
0054 return nullptr;
0055 }
0056
0057
0058 static inline const G4AffineTransform& GetTransform(const G4TouchableHistory* touchable, int depth) {
0059 return touchable->GetHistory()->GetTransform(touchable->GetHistory()->GetDepth() - depth);
0060 }
0061
0062
0063
0064
0065
0066
0067
0068 std::tuple<const G4VPhysicalVolume*, int> GetSensitiveVolume(const G4VTouchable* touchable) {
0069 int depth = touchable->GetHistoryDepth();
0070 for (int level = 0; level < depth; ++level) {
0071 const G4VPhysicalVolume* volume = touchable->GetVolume(level);
0072 if (volume->GetLogicalVolume()->GetSensitiveDetector() != nullptr) {
0073 return std::make_tuple(volume, level);
0074 }
0075 }
0076 return std::tuple<const G4VPhysicalVolume*, int>(nullptr, 0);
0077 }
0078
0079
0080 TrackingMaterialProducer::TrackingMaterialProducer(const edm::ParameterSet& iPSet) {
0081 edm::ParameterSet config = iPSet.getParameter<edm::ParameterSet>("TrackingMaterialProducer");
0082 m_selectedNames = config.getParameter<std::vector<std::string> >("SelectedVolumes");
0083 m_primaryTracks = config.getParameter<bool>("PrimaryTracksOnly");
0084 m_txtOutFile = config.getUntrackedParameter<std::string>("txtOutFile");
0085 m_hgcalzfront = config.getParameter<double>("hgcalzfront");
0086 m_tracks = nullptr;
0087 m_track_volume = nullptr;
0088
0089 produces<std::vector<MaterialAccountingTrack> >();
0090 output_file_ = new TFile("radLen_vs_eta_fromProducer.root", "RECREATE");
0091 output_file_->cd();
0092 radLen_vs_eta_ = new TProfile("radLen", "radLen", 250., -5., 5., 0, 10.);
0093
0094
0095 isHGCal = false;
0096
0097 if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "HGCal") != m_selectedNames.end()) {
0098 isHGCal = true;
0099 }
0100
0101 isHFNose = false;
0102 if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "HFNose") != m_selectedNames.end()) {
0103 isHFNose = true;
0104 }
0105 if (isHGCal or isHFNose) {
0106 outVolumeZpositionTxt.open(m_txtOutFile.c_str(), std::ios::out);
0107 }
0108 }
0109
0110
0111 TrackingMaterialProducer::~TrackingMaterialProducer(void) {}
0112
0113
0114 void TrackingMaterialProducer::update(const EndOfJob* event) {
0115 radLen_vs_eta_->Write();
0116 output_file_->Close();
0117 }
0118
0119 void TrackingMaterialProducer::update(const BeginOfJob* event) {
0120
0121 LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: List of the selected volumes: " << std::endl;
0122 for (std::vector<std::string>::const_iterator volume_name = m_selectedNames.begin();
0123 volume_name != m_selectedNames.end();
0124 ++volume_name) {
0125 const G4LogicalVolume* volume = GetVolume(*volume_name);
0126 if (volume) {
0127 LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: " << *volume_name << std::endl;
0128 m_selectedVolumes.push_back(volume);
0129 } else {
0130
0131 std::cerr << "TrackingMaterialProducer::update(const BeginOfJob*): WARNING: selected volume \"" << *volume_name
0132 << "\" not found in geometry " << std::endl;
0133 }
0134 }
0135 }
0136
0137
0138 void TrackingMaterialProducer::update(const BeginOfEvent* event) {
0139 m_tracks = new std::vector<MaterialAccountingTrack>();
0140 }
0141
0142
0143 void TrackingMaterialProducer::update(const BeginOfTrack* event) {
0144 m_track.reset();
0145 m_track_volume = nullptr;
0146
0147
0148 G4Track* track = const_cast<G4Track*>((*event)());
0149 if (m_primaryTracks and track->GetParentID() != 0) {
0150 track->SetTrackStatus(fStopAndKill);
0151 }
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165 if (isHGCal && track->GetTrackStatus() != fStopAndKill && fabs(track->GetMomentum().eta()) > outerHGCalEta &&
0166 fabs(track->GetMomentum().eta()) < innerHGCalEta) {
0167 if (track->GetMomentum().eta() > 0.) {
0168 outVolumeZpositionTxt << "Air " << m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0 << std::endl;
0169 } else if (track->GetMomentum().eta() <= 0.) {
0170 outVolumeZpositionTxt << "Air " << -m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0 << std::endl;
0171 }
0172 }
0173
0174
0175
0176
0177
0178 if (isHFNose && track->GetTrackStatus() != fStopAndKill && fabs(track->GetMomentum().eta()) > outerHFnoseEta &&
0179 fabs(track->GetMomentum().eta()) < innerHFnoseEta) {
0180 if (track->GetMomentum().eta() > 0.) {
0181 outVolumeZpositionTxt << "Polyethylene " << m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0
0182 << std::endl;
0183 } else if (track->GetMomentum().eta() <= 0.) {
0184 outVolumeZpositionTxt << "Polyethylene " << -m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0
0185 << std::endl;
0186 }
0187 }
0188 }
0189
0190 bool TrackingMaterialProducer::isSelectedFast(const G4TouchableHistory* touchable) {
0191 for (int d = touchable->GetHistoryDepth() - 1; d >= 0; --d) {
0192 if (std::find(m_selectedNames.begin(), m_selectedNames.end(), touchable->GetVolume(d)->GetName()) !=
0193 m_selectedNames.end())
0194 return true;
0195 }
0196 return false;
0197 }
0198
0199
0200 void TrackingMaterialProducer::update(const G4Step* step) {
0201 const G4TouchableHistory* touchable = static_cast<const G4TouchableHistory*>(step->GetTrack()->GetTouchable());
0202 if (not isSelectedFast(touchable)) {
0203 LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer:\t[...] skipping "
0204 << touchable->GetVolume()->GetName() << std::endl;
0205 return;
0206 }
0207
0208
0209 const G4Material* material = touchable->GetVolume()->GetLogicalVolume()->GetMaterial();
0210 double length = step->GetStepLength() / cm;
0211 double X0 = material->GetRadlen() / cm;
0212 double Ne = material->GetElectronDensity() * cm3;
0213 double Xi = Ne / 6.0221415e23 * 0.307075 / 2;
0214 double radiationLengths = length / X0;
0215 double energyLoss = length * Xi / 1000.;
0216
0217
0218 G4ThreeVector globalPosPre = step->GetPreStepPoint()->GetPosition();
0219 G4ThreeVector globalPosPost = step->GetPostStepPoint()->GetPosition();
0220 GlobalPoint globalPositionIn(globalPosPre.x() / cm, globalPosPre.y() / cm, globalPosPre.z() / cm);
0221 GlobalPoint globalPositionOut(globalPosPost.x() / cm, globalPosPost.y() / cm, globalPosPost.z() / cm);
0222
0223 G4StepPoint* prePoint = step->GetPreStepPoint();
0224 G4StepPoint* postPoint = step->GetPostStepPoint();
0225 const CLHEP::Hep3Vector& postPos = postPoint->GetPosition();
0226
0227 if (isHGCal or isHFNose) {
0228
0229
0230
0231 if (postPoint->GetStepStatus() == fGeomBoundary && fabs(postPoint->GetMomentum().eta()) > outerHGCalEta &&
0232 fabs(postPoint->GetMomentum().eta()) < innerHGCalEta) {
0233
0234
0235
0236
0237 outVolumeZpositionTxt << prePoint->GetMaterial()->GetName() << " " << postPos.z() << " "
0238 << postPoint->GetMomentum().eta() << " "
0239 << sqrt(postPos.x() * postPos.x() + postPos.y() * postPos.y()) << " " << 0 << " " << 0
0240 << std::endl;
0241 }
0242
0243 if (postPoint->GetStepStatus() == fGeomBoundary && fabs(postPoint->GetMomentum().eta()) > outerHFnoseEta &&
0244 fabs(postPoint->GetMomentum().eta()) < innerHFnoseEta) {
0245 outVolumeZpositionTxt << prePoint->GetMaterial()->GetName() << " " << postPos.z() << " "
0246 << postPoint->GetMomentum().eta() << " "
0247 << sqrt(postPos.x() * postPos.x() + postPos.y() * postPos.y()) << " " << 0 << " " << 0
0248 << std::endl;
0249 }
0250
0251 }
0252
0253
0254 bool enter_sensitive = false;
0255 bool leave_sensitive = false;
0256 double cosThetaPre = 0.0;
0257 double cosThetaPost = 0.0;
0258 int level = 0;
0259 const G4VPhysicalVolume* sensitive = nullptr;
0260 GlobalPoint position;
0261 std::tie(sensitive, level) = GetSensitiveVolume(touchable);
0262 if (sensitive) {
0263 const G4VSolid& solid = *touchable->GetSolid(level);
0264 const G4AffineTransform& transform = GetTransform(touchable, level);
0265 G4ThreeVector pos = transform.Inverse().TransformPoint(G4ThreeVector(0., 0., 0.));
0266 position = GlobalPoint(pos.x() / cm, pos.y() / cm, pos.z() / cm);
0267
0268 G4ThreeVector localPosPre = transform.TransformPoint(globalPosPre);
0269 EInside statusPre = solid.Inside(localPosPre);
0270 if (statusPre == kSurface) {
0271 enter_sensitive = true;
0272 G4ThreeVector globalDirPre = step->GetPreStepPoint()->GetMomentumDirection();
0273 G4ThreeVector localDirPre = transform.TransformAxis(globalDirPre);
0274 G4ThreeVector normalPre = solid.SurfaceNormal(localPosPre);
0275 cosThetaPre = normalPre.cosTheta(-localDirPre);
0276 }
0277
0278 G4ThreeVector localPosPost = transform.TransformPoint(globalPosPost);
0279 EInside statusPost = solid.Inside(localPosPost);
0280 if (statusPost == kSurface) {
0281 leave_sensitive = true;
0282 G4ThreeVector globalDirPost = step->GetPostStepPoint()->GetMomentumDirection();
0283 G4ThreeVector localDirPost = transform.TransformAxis(globalDirPost);
0284 G4ThreeVector normalPost = solid.SurfaceNormal(localPosPost);
0285 cosThetaPost = normalPost.cosTheta(localDirPost);
0286 }
0287 }
0288
0289
0290 if (enter_sensitive) {
0291 if (m_track_volume != nullptr) {
0292 edm::LogWarning("TrackingMaterialProducer") << "Entering volume " << sensitive << " while inside volume "
0293 << m_track_volume << ". Something is inconsistent";
0294 m_track.reset();
0295 }
0296 m_track_volume = sensitive;
0297 m_track.enterDetector(position, cosThetaPre);
0298 }
0299 m_track.step(MaterialAccountingStep(length, radiationLengths, energyLoss, globalPositionIn, globalPositionOut));
0300 if (leave_sensitive) {
0301 if (m_track_volume != sensitive) {
0302 edm::LogWarning("TrackingMaterialProducer") << "Leaving volume " << sensitive << " while inside volume "
0303 << m_track_volume << ". Something is inconsistent";
0304 m_track.reset();
0305 } else
0306 m_track.leaveDetector(cosThetaPost);
0307 m_track_volume = nullptr;
0308 }
0309 if (sensitive)
0310 LogInfo("TrackingMaterialProducer") << "Track was near sensitive volume " << sensitive->GetName() << std::endl;
0311 else
0312 LogInfo("TrackingMaterialProducer") << "Track was near non-sensitive volume " << touchable->GetVolume()->GetName()
0313 << std::endl;
0314 LogInfo("TrackingMaterialProducer") << "Step length: " << length << " cm\n"
0315 << "globalPreStep(r,z): (" << globalPositionIn.perp() << ", "
0316 << globalPositionIn.z() << ") cm\n"
0317 << "globalPostStep(r,z): (" << globalPositionOut.perp() << ", "
0318 << globalPositionOut.z() << ") cm\n"
0319 << "position(r,z): (" << position.perp() << ", " << position.z() << ") cm\n"
0320 << "Radiation lengths: " << radiationLengths << " \t\t(X0: " << X0
0321 << " cm)\n"
0322 << "Energy loss: " << energyLoss << " MeV \t(Xi: " << Xi
0323 << " MeV/cm)\n"
0324 << "Track was " << (enter_sensitive ? "entering " : "in none ")
0325 << "sensitive volume\n"
0326 << "Track was " << (leave_sensitive ? "leaving " : "in none ")
0327 << "sensitive volume\n";
0328 }
0329
0330
0331 void TrackingMaterialProducer::update(const EndOfTrack* event) {
0332 const G4Track* track = (*event)();
0333 if (m_primaryTracks and track->GetParentID() != 0)
0334 return;
0335
0336 radLen_vs_eta_->Fill(track->GetMomentum().eta(), m_track.summary().radiationLengths());
0337 m_tracks->push_back(m_track);
0338
0339
0340 LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: this track took " << m_track.steps().size()
0341 << " steps, and passed through " << m_track.detectors().size()
0342 << " sensitive detectors" << std::endl;
0343 LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: track length: " << m_track.summary().length()
0344 << " cm" << std::endl;
0345 LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: radiation lengths: "
0346 << m_track.summary().radiationLengths() << std::endl;
0347 LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: energy loss: "
0348 << m_track.summary().energyLoss() << " MeV" << std::endl;
0349 }
0350
0351
0352 void TrackingMaterialProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0353
0354 std::unique_ptr<std::vector<MaterialAccountingTrack> > tracks(m_tracks);
0355 iEvent.put(std::move(tracks));
0356 m_tracks = nullptr;
0357 }
0358
0359
0360 bool TrackingMaterialProducer::isSelected(const G4VTouchable* touchable) {
0361 for (size_t i = 0; i < m_selectedVolumes.size(); ++i)
0362 if (m_selectedVolumes[i]->IsAncestor(touchable->GetVolume()) or
0363 m_selectedVolumes[i] == touchable->GetVolume()->GetLogicalVolume())
0364 return true;
0365
0366 return false;
0367 }
0368
0369
0370
0371 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0372 #include "FWCore/Framework/interface/MakerMacros.h"
0373 DEFINE_SIMWATCHER(TrackingMaterialProducer);