Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:04 

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 <DD4hep/Filter.h>
0031 
0032 #include "TrackingMaterialProducer.h"
0033 
0034 // Uncomment the following #define directive to have the full list of
0035 // volumes known to G4 printed to LogInfo("TrackingMaterialProducer")
0036 
0037 #define DEBUG_G4_VOLUMES
0038 
0039 using namespace CLHEP;
0040 using edm::LogInfo;
0041 
0042 // missing from GEANT4 < 9.0 : G4LogicalVolumeStore::GetVolume( name )
0043 static const G4LogicalVolume* GetVolume(const std::string& name) {
0044   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0045 
0046 #ifdef DEBUG_G4_VOLUMES
0047   for (G4LogicalVolumeStore::const_iterator volume = lvs->begin(); volume != lvs->end(); ++volume)
0048     LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: G4 registered volumes " << (*volume)->GetName()
0049                                         << std::endl;
0050 #endif
0051 
0052   for (G4LogicalVolumeStore::const_iterator volume = lvs->begin(); volume != lvs->end(); ++volume) {
0053     if ((const std::string)(dd4hep::dd::noNamespace((*volume)->GetName())) == name)
0054       return (*volume);
0055   }
0056   return nullptr;
0057 }
0058 
0059 // missing from GEANT4 : G4TouchableHistory::GetTransform( depth )
0060 static inline const G4AffineTransform& GetTransform(const G4TouchableHistory* touchable, int depth) {
0061   return touchable->GetHistory()->GetTransform(touchable->GetHistory()->GetDepth() - depth);
0062 }
0063 
0064 // navigate up the hierarchy of volumes until one with an attached sensitive detector is found
0065 // return a tuple holding
0066 //   - pointer to the first (deepest) sensitive G4VPhysicalVolume
0067 //   - how may steps up in the hierarchy it is (0 is the starting volume)
0068 // if no sensitive detector is found, return a NULL pointer and 0
0069 
0070 std::tuple<const G4VPhysicalVolume*, int> GetSensitiveVolume(const G4VTouchable* touchable) {
0071   int depth = touchable->GetHistoryDepth();
0072   for (int level = 0; level < depth; ++level) {  // 0 is self
0073     const G4VPhysicalVolume* volume = touchable->GetVolume(level);
0074     if (volume->GetLogicalVolume()->GetSensitiveDetector() != nullptr) {
0075       return std::make_tuple(volume, level);
0076     }
0077   }
0078   return std::tuple<const G4VPhysicalVolume*, int>(nullptr, 0);
0079 }
0080 
0081 //-------------------------------------------------------------------------
0082 TrackingMaterialProducer::TrackingMaterialProducer(const edm::ParameterSet& iPSet) {
0083   edm::ParameterSet config = iPSet.getParameter<edm::ParameterSet>("TrackingMaterialProducer");
0084   m_selectedNames = config.getParameter<std::vector<std::string> >("SelectedVolumes");
0085   m_primaryTracks = config.getParameter<bool>("PrimaryTracksOnly");
0086   m_txtOutFile = config.getUntrackedParameter<std::string>("txtOutFile");
0087   m_hgcalzfront = config.getParameter<double>("hgcalzfront");
0088   m_tracks = nullptr;
0089   m_track_volume = nullptr;
0090 
0091   produces<std::vector<MaterialAccountingTrack> >();
0092   output_file_ = new TFile("radLen_vs_eta_fromProducer.root", "RECREATE");
0093   output_file_->cd();
0094   radLen_vs_eta_ = new TProfile("radLen", "radLen", 250., -5., 5., 0, 10.);
0095 
0096   //Check if HGCal volumes are selected
0097   isHGCal = false;
0098   //if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "CALOECTSRear") != m_selectedNames.end()) {
0099   if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "HGCal") != m_selectedNames.end()) {
0100     isHGCal = true;
0101   }
0102   //Check if HFNose volumes are selected
0103   isHFNose = false;
0104   if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "HFNose") != m_selectedNames.end()) {
0105     isHFNose = true;
0106   }
0107   if (isHGCal or isHFNose) {
0108     outVolumeZpositionTxt.open(m_txtOutFile.c_str(), std::ios::out);
0109   }
0110 }
0111 
0112 //-------------------------------------------------------------------------
0113 TrackingMaterialProducer::~TrackingMaterialProducer(void) {}
0114 
0115 //-------------------------------------------------------------------------
0116 void TrackingMaterialProducer::update(const EndOfJob* event) {
0117   radLen_vs_eta_->Write();
0118   output_file_->Close();
0119 }
0120 //-------------------------------------------------------------------------
0121 void TrackingMaterialProducer::update(const BeginOfJob* event) {
0122   // INFO
0123   LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: List of the selected volumes: " << std::endl;
0124   for (std::vector<std::string>::const_iterator volume_name = m_selectedNames.begin();
0125        volume_name != m_selectedNames.end();
0126        ++volume_name) {
0127     const G4LogicalVolume* volume = GetVolume(*volume_name);
0128     if (volume) {
0129       LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: " << *volume_name << std::endl;
0130       m_selectedVolumes.push_back(volume);
0131     } else {
0132       // FIXME: throw an exception ?
0133       std::cerr << "TrackingMaterialProducer::update(const BeginOfJob*): WARNING: selected volume \"" << *volume_name
0134                 << "\" not found in geometry " << std::endl;
0135     }
0136   }
0137 }
0138 
0139 //-------------------------------------------------------------------------
0140 void TrackingMaterialProducer::update(const BeginOfEvent* event) {
0141   m_tracks = new std::vector<MaterialAccountingTrack>();
0142 }
0143 
0144 //-------------------------------------------------------------------------
0145 void TrackingMaterialProducer::update(const BeginOfTrack* event) {
0146   m_track.reset();
0147   m_track_volume = nullptr;
0148 
0149   // prevent secondary tracks from propagating
0150   G4Track* track = const_cast<G4Track*>((*event)());
0151   if (m_primaryTracks and track->GetParentID() != 0) {
0152     track->SetTrackStatus(fStopAndKill);
0153   }
0154 
0155   //For the HGCal case:
0156   //In the beginning of each track, the track will first hit an HGCAL volume and it will
0157   //save the upper z volume boundary. So, the low boundary of the first
0158   //volume is never saved. Here we give the low boundary of the first volume.
0159   //This can be found by asking first to run not on 'HGCal' volume below but
0160   //on 'CALOECTSRear', which at the moment of this writing it contains
0161   //HGCalService, HGCal and thermal screen. You should run Fireworks to
0162   //check if these naming conventions and volumes are valid in the future.
0163   //Then, check the VolumesZPosition.txt file to see where CEService ends and
0164   //put that number in hgcalzfront. Keep in mind to run on the desired volume above here:
0165   //https://github.com/cms-sw/cmssw/blob/master/SimTracker/TrackerMaterialAnalysis/plugins/TrackingMaterialProducer.cc#L95
0166   //and to replace the volume name of the material first hit at the file creation line below
0167   if (isHGCal && track->GetTrackStatus() != fStopAndKill && fabs(track->GetMomentum().eta()) > outerHGCalEta &&
0168       fabs(track->GetMomentum().eta()) < innerHGCalEta) {
0169     if (track->GetMomentum().eta() > 0.) {
0170       outVolumeZpositionTxt << "Air " << m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0 << std::endl;
0171     } else if (track->GetMomentum().eta() <= 0.) {
0172       outVolumeZpositionTxt << "Air " << -m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0 << std::endl;
0173     }
0174   }
0175 
0176   //For the HFnose case:
0177   //restrict the outher radius to eta 3.3 since there is HGCAL shadowing
0178   //restrict the innner radius to eta 4 since it's non projective
0179 
0180   if (isHFNose && track->GetTrackStatus() != fStopAndKill && fabs(track->GetMomentum().eta()) > outerHFnoseEta &&
0181       fabs(track->GetMomentum().eta()) < innerHFnoseEta) {
0182     if (track->GetMomentum().eta() > 0.) {
0183       outVolumeZpositionTxt << "Polyethylene " << m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0
0184                             << std::endl;
0185     } else if (track->GetMomentum().eta() <= 0.) {
0186       outVolumeZpositionTxt << "Polyethylene " << -m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0
0187                             << std::endl;
0188     }
0189   }
0190 }
0191 
0192 bool TrackingMaterialProducer::isSelectedFast(const G4TouchableHistory* touchable) {
0193   for (int d = touchable->GetHistoryDepth() - 1; d >= 0; --d) {
0194     if (std::find(m_selectedNames.begin(),
0195                   m_selectedNames.end(),
0196                   (std::string)(dd4hep::dd::noNamespace(touchable->GetVolume(d)->GetName()))) != m_selectedNames.end())
0197       return true;
0198   }
0199   return false;
0200 }
0201 
0202 //-------------------------------------------------------------------------
0203 void TrackingMaterialProducer::update(const G4Step* step) {
0204   const G4TouchableHistory* touchable = static_cast<const G4TouchableHistory*>(step->GetTrack()->GetTouchable());
0205   if (not isSelectedFast(touchable)) {
0206     LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer:\t[...] skipping "
0207                                         << touchable->GetVolume()->GetName() << std::endl;
0208     return;
0209   }
0210 
0211   // material and step proterties
0212   const G4Material* material = touchable->GetVolume()->GetLogicalVolume()->GetMaterial();
0213   double length = step->GetStepLength() / cm;        // mm -> cm
0214   double X0 = material->GetRadlen() / cm;            // mm -> cm
0215   double Ne = material->GetElectronDensity() * cm3;  // 1/mm3 -> 1/cm3
0216   double Xi = Ne / 6.0221415e23 * 0.307075 / 2;      // MeV / cm
0217   double radiationLengths = length / X0;             //
0218   double energyLoss = length * Xi / 1000.;           // GeV
0219   //double energyLoss = step->GetDeltaEnergy()/MeV;  should we use this??
0220 
0221   G4ThreeVector globalPosPre = step->GetPreStepPoint()->GetPosition();
0222   G4ThreeVector globalPosPost = step->GetPostStepPoint()->GetPosition();
0223   GlobalPoint globalPositionIn(globalPosPre.x() / cm, globalPosPre.y() / cm, globalPosPre.z() / cm);      // mm -> cm
0224   GlobalPoint globalPositionOut(globalPosPost.x() / cm, globalPosPost.y() / cm, globalPosPost.z() / cm);  // mm -> cm
0225 
0226   G4StepPoint* prePoint = step->GetPreStepPoint();
0227   G4StepPoint* postPoint = step->GetPostStepPoint();
0228   const CLHEP::Hep3Vector& postPos = postPoint->GetPosition();
0229   //Go below only in HGCal case
0230   if (isHGCal or isHFNose) {
0231     //A step never spans across boundaries: geometry or physics define the end points
0232     //If the step is limited by a boundary, the post-step point stands on the
0233     //boundary and it logically belongs to the next volume.
0234     if (postPoint->GetStepStatus() == fGeomBoundary && fabs(postPoint->GetMomentum().eta()) > outerHGCalEta &&
0235         fabs(postPoint->GetMomentum().eta()) < innerHGCalEta) {
0236       //Post point position is the low z edge of the new volume, or the upper for the prepoint volume.
0237       //So, premat - postz - posteta - postR - premattotalenergylossEtable - premattotalenergylossEfull
0238       //Observe the two zeros at the end which in the past where set to emCalculator.GetDEDX and
0239       //emCalculator.ComputeTotalDEDX but decided not to be used. Will save the structure though for the script.
0240       outVolumeZpositionTxt << prePoint->GetMaterial()->GetName() << " " << postPos.z() << " "
0241                             << postPoint->GetMomentum().eta() << " "
0242                             << sqrt(postPos.x() * postPos.x() + postPos.y() * postPos.y()) << " " << 0 << " " << 0
0243                             << std::endl;
0244     }
0245 
0246     if (postPoint->GetStepStatus() == fGeomBoundary && fabs(postPoint->GetMomentum().eta()) > outerHFnoseEta &&
0247         fabs(postPoint->GetMomentum().eta()) < innerHFnoseEta) {
0248       outVolumeZpositionTxt << prePoint->GetMaterial()->GetName() << " " << postPos.z() << " "
0249                             << postPoint->GetMomentum().eta() << " "
0250                             << sqrt(postPos.x() * postPos.x() + postPos.y() * postPos.y()) << " " << 0 << " " << 0
0251                             << std::endl;
0252     }
0253 
0254   }  //end of isHGCal  or HFnose if
0255 
0256   // check for a sensitive detector
0257   bool enter_sensitive = false;
0258   bool leave_sensitive = false;
0259   double cosThetaPre = 0.0;
0260   double cosThetaPost = 0.0;
0261   int level = 0;
0262   const G4VPhysicalVolume* sensitive = nullptr;
0263   GlobalPoint position;
0264   std::tie(sensitive, level) = GetSensitiveVolume(touchable);
0265   if (sensitive) {
0266     const G4VSolid& solid = *touchable->GetSolid(level);
0267     const G4AffineTransform& transform = GetTransform(touchable, level);
0268     G4ThreeVector pos = transform.Inverse().TransformPoint(G4ThreeVector(0., 0., 0.));
0269     position = GlobalPoint(pos.x() / cm, pos.y() / cm, pos.z() / cm);  // mm -> cm
0270 
0271     G4ThreeVector localPosPre = transform.TransformPoint(globalPosPre);
0272     EInside statusPre = solid.Inside(localPosPre);
0273     if (statusPre == kSurface) {
0274       enter_sensitive = true;
0275       G4ThreeVector globalDirPre = step->GetPreStepPoint()->GetMomentumDirection();
0276       G4ThreeVector localDirPre = transform.TransformAxis(globalDirPre);
0277       G4ThreeVector normalPre = solid.SurfaceNormal(localPosPre);
0278       cosThetaPre = normalPre.cosTheta(-localDirPre);
0279     }
0280 
0281     G4ThreeVector localPosPost = transform.TransformPoint(globalPosPost);
0282     EInside statusPost = solid.Inside(localPosPost);
0283     if (statusPost == kSurface) {
0284       leave_sensitive = true;
0285       G4ThreeVector globalDirPost = step->GetPostStepPoint()->GetMomentumDirection();
0286       G4ThreeVector localDirPost = transform.TransformAxis(globalDirPost);
0287       G4ThreeVector normalPost = solid.SurfaceNormal(localPosPost);
0288       cosThetaPost = normalPost.cosTheta(localDirPost);
0289     }
0290   }
0291 
0292   // update track accounting
0293   if (enter_sensitive) {
0294     if (m_track_volume != nullptr) {
0295       edm::LogWarning("TrackingMaterialProducer") << "Entering volume " << sensitive << " while inside volume "
0296                                                   << m_track_volume << ". Something is inconsistent";
0297       m_track.reset();
0298     }
0299     m_track_volume = sensitive;
0300     m_track.enterDetector(position, cosThetaPre);
0301   }
0302   m_track.step(MaterialAccountingStep(length, radiationLengths, energyLoss, globalPositionIn, globalPositionOut));
0303   if (leave_sensitive) {
0304     if (m_track_volume != sensitive) {
0305       edm::LogWarning("TrackingMaterialProducer") << "Leaving volume " << sensitive << " while inside volume "
0306                                                   << m_track_volume << ". Something is inconsistent";
0307       m_track.reset();
0308     } else
0309       m_track.leaveDetector(cosThetaPost);
0310     m_track_volume = nullptr;
0311   }
0312   if (sensitive)
0313     LogInfo("TrackingMaterialProducer") << "Track was near sensitive     volume " << sensitive->GetName() << std::endl;
0314   else
0315     LogInfo("TrackingMaterialProducer") << "Track was near non-sensitive volume " << touchable->GetVolume()->GetName()
0316                                         << std::endl;
0317   LogInfo("TrackingMaterialProducer") << "Step length:             " << length << " cm\n"
0318                                       << "globalPreStep(r,z): (" << globalPositionIn.perp() << ", "
0319                                       << globalPositionIn.z() << ") cm\n"
0320                                       << "globalPostStep(r,z): (" << globalPositionOut.perp() << ", "
0321                                       << globalPositionOut.z() << ") cm\n"
0322                                       << "position(r,z): (" << position.perp() << ", " << position.z() << ") cm\n"
0323                                       << "Radiation lengths:       " << radiationLengths << " \t\t(X0: " << X0
0324                                       << " cm)\n"
0325                                       << "Energy loss:             " << energyLoss << " MeV  \t(Xi: " << Xi
0326                                       << " MeV/cm)\n"
0327                                       << "Track was " << (enter_sensitive ? "entering " : "in none ")
0328                                       << "sensitive volume\n"
0329                                       << "Track was " << (leave_sensitive ? "leaving  " : "in none ")
0330                                       << "sensitive volume\n";
0331 }
0332 
0333 //-------------------------------------------------------------------------
0334 void TrackingMaterialProducer::update(const EndOfTrack* event) {
0335   const G4Track* track = (*event)();
0336   if (m_primaryTracks and track->GetParentID() != 0)
0337     return;
0338 
0339   radLen_vs_eta_->Fill(track->GetMomentum().eta(), m_track.summary().radiationLengths());
0340   m_tracks->push_back(m_track);
0341 
0342   // LogInfo
0343   LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: this track took " << m_track.steps().size()
0344                                       << " steps, and passed through " << m_track.detectors().size()
0345                                       << " sensitive detectors" << std::endl;
0346   LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: track length:       " << m_track.summary().length()
0347                                       << " cm" << std::endl;
0348   LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: radiation lengths: "
0349                                       << m_track.summary().radiationLengths() << std::endl;
0350   LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: energy loss:        "
0351                                       << m_track.summary().energyLoss() << " MeV" << std::endl;
0352 }
0353 
0354 //-------------------------------------------------------------------------
0355 void TrackingMaterialProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0356   // transfer ownership to the Event
0357   std::unique_ptr<std::vector<MaterialAccountingTrack> > tracks(m_tracks);
0358   iEvent.put(std::move(tracks));
0359   m_tracks = nullptr;
0360 }
0361 
0362 //-------------------------------------------------------------------------
0363 bool TrackingMaterialProducer::isSelected(const G4VTouchable* touchable) {
0364   for (size_t i = 0; i < m_selectedVolumes.size(); ++i)
0365     if (m_selectedVolumes[i]->IsAncestor(touchable->GetVolume()) or
0366         m_selectedVolumes[i] == touchable->GetVolume()->GetLogicalVolume())
0367       return true;
0368 
0369   return false;
0370 }
0371 
0372 //-------------------------------------------------------------------------
0373 // define as a plugin
0374 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0375 #include "FWCore/Framework/interface/MakerMacros.h"
0376 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!