Back to home page

Project CMSSW displayed by LXR
 
 

    

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 // GEANT4
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 // Uncomment the following #define directive to have the full list of
0033 // volumes known to G4 printed to LogInfo("TrackingMaterialProducer")
0034 
0035 #define DEBUG_G4_VOLUMES
0036 
0037 using namespace CLHEP;
0038 using edm::LogInfo;
0039 
0040 // missing from GEANT4 < 9.0 : G4LogicalVolumeStore::GetVolume( name )
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 // missing from GEANT4 : G4TouchableHistory::GetTransform( depth )
0058 static inline const G4AffineTransform& GetTransform(const G4TouchableHistory* touchable, int depth) {
0059   return touchable->GetHistory()->GetTransform(touchable->GetHistory()->GetDepth() - depth);
0060 }
0061 
0062 // navigate up the hierarchy of volumes until one with an attached sensitive detector is found
0063 // return a tuple holding
0064 //   - pointer to the first (deepest) sensitive G4VPhysicalVolume
0065 //   - how may steps up in the hierarchy it is (0 is the starting volume)
0066 // if no sensitive detector is found, return a NULL pointer and 0
0067 
0068 std::tuple<const G4VPhysicalVolume*, int> GetSensitiveVolume(const G4VTouchable* touchable) {
0069   int depth = touchable->GetHistoryDepth();
0070   for (int level = 0; level < depth; ++level) {  // 0 is self
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   //Check if HGCal volumes are selected
0095   isHGCal = false;
0096   //if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "CALOECTSRear") != m_selectedNames.end()) {
0097   if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "HGCal") != m_selectedNames.end()) {
0098     isHGCal = true;
0099   }
0100   //Check if HFNose volumes are selected
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   // INFO
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       // FIXME: throw an exception ?
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   // prevent secondary tracks from propagating
0148   G4Track* track = const_cast<G4Track*>((*event)());
0149   if (m_primaryTracks and track->GetParentID() != 0) {
0150     track->SetTrackStatus(fStopAndKill);
0151   }
0152 
0153   //For the HGCal case:
0154   //In the beginning of each track, the track will first hit an HGCAL volume and it will
0155   //save the upper z volume boundary. So, the low boundary of the first
0156   //volume is never saved. Here we give the low boundary of the first volume.
0157   //This can be found by asking first to run not on 'HGCal' volume below but
0158   //on 'CALOECTSRear', which at the moment of this writing it contains
0159   //HGCalService, HGCal and thermal screen. You should run Fireworks to
0160   //check if these naming conventions and volumes are valid in the future.
0161   //Then, check the VolumesZPosition.txt file to see where CEService ends and
0162   //put that number in hgcalzfront. Keep in mind to run on the desired volume above here:
0163   //https://github.com/cms-sw/cmssw/blob/master/SimTracker/TrackerMaterialAnalysis/plugins/TrackingMaterialProducer.cc#L95
0164   //and to replace the volume name of the material first hit at the file creation line below
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   //For the HFnose case:
0175   //restrict the outher radius to eta 3.3 since there is HGCAL shadowing
0176   //restrict the innner radius to eta 4 since it's non projective
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   // material and step proterties
0209   const G4Material* material = touchable->GetVolume()->GetLogicalVolume()->GetMaterial();
0210   double length = step->GetStepLength() / cm;        // mm -> cm
0211   double X0 = material->GetRadlen() / cm;            // mm -> cm
0212   double Ne = material->GetElectronDensity() * cm3;  // 1/mm3 -> 1/cm3
0213   double Xi = Ne / 6.0221415e23 * 0.307075 / 2;      // MeV / cm
0214   double radiationLengths = length / X0;             //
0215   double energyLoss = length * Xi / 1000.;           // GeV
0216   //double energyLoss = step->GetDeltaEnergy()/MeV;  should we use this??
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);      // mm -> cm
0221   GlobalPoint globalPositionOut(globalPosPost.x() / cm, globalPosPost.y() / cm, globalPosPost.z() / cm);  // mm -> cm
0222 
0223   G4StepPoint* prePoint = step->GetPreStepPoint();
0224   G4StepPoint* postPoint = step->GetPostStepPoint();
0225   const CLHEP::Hep3Vector& postPos = postPoint->GetPosition();
0226   //Go below only in HGCal case
0227   if (isHGCal or isHFNose) {
0228     //A step never spans across boundaries: geometry or physics define the end points
0229     //If the step is limited by a boundary, the post-step point stands on the
0230     //boundary and it logically belongs to the next volume.
0231     if (postPoint->GetStepStatus() == fGeomBoundary && fabs(postPoint->GetMomentum().eta()) > outerHGCalEta &&
0232         fabs(postPoint->GetMomentum().eta()) < innerHGCalEta) {
0233       //Post point position is the low z edge of the new volume, or the upper for the prepoint volume.
0234       //So, premat - postz - posteta - postR - premattotalenergylossEtable - premattotalenergylossEfull
0235       //Observe the two zeros at the end which in the past where set to emCalculator.GetDEDX and
0236       //emCalculator.ComputeTotalDEDX but decided not to be used. Will save the structure though for the script.
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   }  //end of isHGCal  or HFnose if
0252 
0253   // check for a sensitive detector
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);  // mm -> 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   // update track accounting
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   // LogInfo
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   // transfer ownership to the Event
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 // define as a plugin
0371 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0372 #include "FWCore/Framework/interface/MakerMacros.h"
0373 DEFINE_SIMWATCHER(TrackingMaterialProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!