TrackingMaterialProducer Class Reference

#include <TrackingMaterialProducer.h>

Inheritance diagram for TrackingMaterialProducer:

Public Member Functions
	TrackingMaterialProducer (const edm::ParameterSet &)
	~TrackingMaterialProducer () override
Public Member Functions inherited from SimProducer
const SimProducer &	operator= (const SimProducer &)=delete
void	registerProducts (edm::ProducesCollector producesCollector)
	SimProducer ()
	SimProducer (const SimProducer &)=delete
Public Member Functions inherited from SimWatcher
virtual void	beginRun (edm::EventSetup const &)
bool	isMT () const
const SimWatcher &	operator= (const SimWatcher &)=delete
virtual void	registerConsumes (edm::ConsumesCollector)
	SimWatcher ()
	SimWatcher (const SimWatcher &)=delete
virtual	~SimWatcher ()
Public Member Functions inherited from Observer< const BeginOfJob * >
	Observer ()
void	slotForUpdate (const BeginOfJob *iT)
virtual	~Observer ()
Public Member Functions inherited from Observer< const EndOfJob * >
	Observer ()
void	slotForUpdate (const EndOfJob *iT)
virtual	~Observer ()
Public Member Functions inherited from Observer< const BeginOfEvent * >
	Observer ()
void	slotForUpdate (const BeginOfEvent *iT)
virtual	~Observer ()
Public Member Functions inherited from Observer< const BeginOfTrack * >
	Observer ()
void	slotForUpdate (const BeginOfTrack *iT)
virtual	~Observer ()
Public Member Functions inherited from Observer< const G4Step * >
	Observer ()
void	slotForUpdate (const G4Step *iT)
virtual	~Observer ()
Public Member Functions inherited from Observer< const EndOfTrack * >
	Observer ()
void	slotForUpdate (const EndOfTrack *iT)
virtual	~Observer ()

Private Member Functions
bool	isSelected (const G4VTouchable *touch)
bool	isSelectedFast (const G4TouchableHistory *touch)
void	produce (edm::Event &, const edm::EventSetup &) override
void	update (const BeginOfJob *) override
	This routine will be called when the appropriate signal arrives. More...
void	update (const BeginOfEvent *) override
	This routine will be called when the appropriate signal arrives. More...
void	update (const BeginOfTrack *) override
	This routine will be called when the appropriate signal arrives. More...
void	update (const G4Step *) override
	This routine will be called when the appropriate signal arrives. More...
void	update (const EndOfTrack *) override
	This routine will be called when the appropriate signal arrives. More...
void	update (const EndOfJob *) override
	This routine will be called when the appropriate signal arrives. More...

Private Attributes
bool	isHFNose
bool	isHGCal
double	m_hgcalzfront
bool	m_primaryTracks
std::vector< std::string >	m_selectedNames
std::vector< const G4LogicalVolume * >	m_selectedVolumes
MaterialAccountingTrack	m_track
const G4VPhysicalVolume *	m_track_volume
std::vector < MaterialAccountingTrack > *	m_tracks
std::string	m_txtOutFile
TFile *	output_file_
std::ofstream	outVolumeZpositionTxt
TProfile *	radLen_vs_eta_

Static Private Attributes
static constexpr float	innerHFnoseEta = 4.
static constexpr float	innerHGCalEta = 2.4
static constexpr float	outerHFnoseEta = 3.3
static constexpr float	outerHGCalEta = 2.0

Additional Inherited Members
Protected Member Functions inherited from SimProducer
template<class T >
void	produces ()
template<class T >
void	produces (const std::string &instanceName)
Protected Member Functions inherited from SimWatcher
void	setMT (bool val)

Detailed Description

Definition at line 33 of file TrackingMaterialProducer.h.

Constructor & Destructor Documentation

TrackingMaterialProducer::TrackingMaterialProducer

(

const edm::ParameterSet &

iPSet

)

Definition at line 80 of file TrackingMaterialProducer.cc.

References submitPVResolutionJobs::config, spr::find(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), or, submitPVResolutionJobs::out, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                               {
  edm::ParameterSet config = iPSet.getParameter<edm::ParameterSet>("TrackingMaterialProducer");
  m_selectedNames = config.getParameter<std::vector<std::string> >("SelectedVolumes");
  m_primaryTracks = config.getParameter<bool>("PrimaryTracksOnly");
  m_txtOutFile = config.getUntrackedParameter<std::string>("txtOutFile");
  m_hgcalzfront = config.getParameter<double>("hgcalzfront");
  m_tracks = nullptr;
  m_track_volume = nullptr;

  produces<std::vector<MaterialAccountingTrack> >();
  output_file_ = new TFile("radLen_vs_eta_fromProducer.root", "RECREATE");
  output_file_->cd();
  radLen_vs_eta_ = new TProfile("radLen", "radLen", 250., -5., 5., 0, 10.);

  //Check if HGCal volumes are selected
  isHGCal = false;
  //if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "CALOECTSRear") != m_selectedNames.end()) {
  if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "HGCal") != m_selectedNames.end()) {
    isHGCal = true;
  }
  //Check if HFNose volumes are selected
  isHFNose = false;
  if (std::find(m_selectedNames.begin(), m_selectedNames.end(), "HFNose") != m_selectedNames.end()) {
    isHFNose = true;
  }
  if (isHGCal or isHFNose) {
    outVolumeZpositionTxt.open(m_txtOutFile.c_str(), std::ios::out);
  }
}

TrackingMaterialProducer::~TrackingMaterialProducer ( void  )

override

Definition at line 111 of file TrackingMaterialProducer.cc.

{}

Member Function Documentation

bool TrackingMaterialProducer::isSelected ( const G4VTouchable *  touch )

private

Definition at line 360 of file TrackingMaterialProducer.cc.

References mps_fire::i, and or.

                                                                       {
  for (size_t i = 0; i < m_selectedVolumes.size(); ++i)
    if (m_selectedVolumes[i]->IsAncestor(touchable->GetVolume()) or
        m_selectedVolumes[i] == touchable->GetVolume()->GetLogicalVolume())
      return true;

  return false;
}

bool TrackingMaterialProducer::isSelectedFast ( const G4TouchableHistory *  touch )

private

Definition at line 190 of file TrackingMaterialProducer.cc.

References ztail::d, and spr::find().

                                                                                 {
  for (int d = touchable->GetHistoryDepth() - 1; d >= 0; --d) {
    if (std::find(m_selectedNames.begin(), m_selectedNames.end(), touchable->GetVolume(d)->GetName()) !=
        m_selectedNames.end())
      return true;
  }
  return false;
}

void TrackingMaterialProducer::produce ( edm::Event &  iEvent, const edm::EventSetup &  iSetup  )

overrideprivatevirtual

Implements SimProducer.

Definition at line 352 of file TrackingMaterialProducer.cc.

References eostools::move(), edm::Event::put(), and tracks.

                                                                                    {
  // transfer ownership to the Event
  std::unique_ptr<std::vector<MaterialAccountingTrack> > tracks(m_tracks);
  iEvent.put(std::move(tracks));
  m_tracks = nullptr;
}

void TrackingMaterialProducer::update ( const BeginOfJob *  )

overrideprivatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfJob * >.

Definition at line 119 of file TrackingMaterialProducer.cc.

References EcnaPython_AdcPeg12_S1_10_R170298_1_0_150_Dee0::cerr, and GetVolume().

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                             {
  // INFO
  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: List of the selected volumes: " << std::endl;
  for (std::vector<std::string>::const_iterator volume_name = m_selectedNames.begin();
       volume_name != m_selectedNames.end();
       ++volume_name) {
    const G4LogicalVolume* volume = GetVolume(*volume_name);
    if (volume) {
      LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: " << *volume_name << std::endl;
      m_selectedVolumes.push_back(volume);
    } else {
      // FIXME: throw an exception ?
      std::cerr << "TrackingMaterialProducer::update(const BeginOfJob*): WARNING: selected volume \"" << *volume_name
                << "\" not found in geometry " << std::endl;
    }
  }
}

void TrackingMaterialProducer::update ( const BeginOfEvent *  )

overrideprivatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfEvent * >.

Definition at line 138 of file TrackingMaterialProducer.cc.

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                               {
  m_tracks = new std::vector<MaterialAccountingTrack>();
}

void TrackingMaterialProducer::update ( const BeginOfTrack *  )

overrideprivatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfTrack * >.

Definition at line 143 of file TrackingMaterialProducer.cc.

References HLT_FULL_cff::track.

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                               {
  m_track.reset();
  m_track_volume = nullptr;

  // prevent secondary tracks from propagating
  G4Track* track = const_cast<G4Track*>((*event)());
  if (m_primaryTracks and track->GetParentID() != 0) {
    track->SetTrackStatus(fStopAndKill);
  }

  //For the HGCal case:
  //In the beginning of each track, the track will first hit an HGCAL volume and it will
  //save the upper z volume boundary. So, the low boundary of the first
  //volume is never saved. Here we give the low boundary of the first volume.
  //This can be found by asking first to run not on 'HGCal' volume below but
  //on 'CALOECTSRear', which at the moment of this writing it contains
  //HGCalService, HGCal and thermal screen. You should run Fireworks to
  //check if these naming conventions and volumes are valid in the future.
  //Then, check the VolumesZPosition.txt file to see where CEService ends and
  //put that number in hgcalzfront. Keep in mind to run on the desired volume above here:
  //https://github.com/cms-sw/cmssw/blob/master/SimTracker/TrackerMaterialAnalysis/plugins/TrackingMaterialProducer.cc#L95
  //and to replace the volume name of the material first hit at the file creation line below
  if (isHGCal && track->GetTrackStatus() != fStopAndKill && fabs(track->GetMomentum().eta()) > outerHGCalEta &&
      fabs(track->GetMomentum().eta()) < innerHGCalEta) {
    if (track->GetMomentum().eta() > 0.) {
      outVolumeZpositionTxt << "Air " << m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0 << std::endl;
    } else if (track->GetMomentum().eta() <= 0.) {
      outVolumeZpositionTxt << "Air " << -m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0 << std::endl;
    }
  }

  //For the HFnose case:
  //restrict the outher radius to eta 3.3 since there is HGCAL shadowing
  //restrict the innner radius to eta 4 since it's non projective

  if (isHFNose && track->GetTrackStatus() != fStopAndKill && fabs(track->GetMomentum().eta()) > outerHFnoseEta &&
      fabs(track->GetMomentum().eta()) < innerHFnoseEta) {
    if (track->GetMomentum().eta() > 0.) {
      outVolumeZpositionTxt << "Polyethylene " << m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0
                            << std::endl;
    } else if (track->GetMomentum().eta() <= 0.) {
      outVolumeZpositionTxt << "Polyethylene " << -m_hgcalzfront << " " << 0 << " " << 0 << " " << 0 << " " << 0
                            << std::endl;
    }
  }
}

void TrackingMaterialProducer::update ( const G4Step *  )

overrideprivatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const G4Step * >.

Definition at line 200 of file TrackingMaterialProducer.cc.

References GetSensitiveVolume(), GetTransform(), testEve_cfg::level, or, PV3DBase< T, PVType, FrameType >::perp(), position, mathSSE::sqrt(), HcalDetIdTransform::transform(), X0, and PV3DBase< T, PVType, FrameType >::z().

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                        {
  const G4TouchableHistory* touchable = static_cast<const G4TouchableHistory*>(step->GetTrack()->GetTouchable());
  if (not isSelectedFast(touchable)) {
    LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer:\t[...] skipping "
                                        << touchable->GetVolume()->GetName() << std::endl;
    return;
  }

  // material and step proterties
  const G4Material* material = touchable->GetVolume()->GetLogicalVolume()->GetMaterial();
  double length = step->GetStepLength() / cm;        // mm -> cm
  double X0 = material->GetRadlen() / cm;            // mm -> cm
  double Ne = material->GetElectronDensity() * cm3;  // 1/mm3 -> 1/cm3
  double Xi = Ne / 6.0221415e23 * 0.307075 / 2;      // MeV / cm
  double radiationLengths = length / X0;             //
  double energyLoss = length * Xi / 1000.;           // GeV
  //double energyLoss = step->GetDeltaEnergy()/MeV;  should we use this??

  G4ThreeVector globalPosPre = step->GetPreStepPoint()->GetPosition();
  G4ThreeVector globalPosPost = step->GetPostStepPoint()->GetPosition();
  GlobalPoint globalPositionIn(globalPosPre.x() / cm, globalPosPre.y() / cm, globalPosPre.z() / cm);      // mm -> cm
  GlobalPoint globalPositionOut(globalPosPost.x() / cm, globalPosPost.y() / cm, globalPosPost.z() / cm);  // mm -> cm

  G4StepPoint* prePoint = step->GetPreStepPoint();
  G4StepPoint* postPoint = step->GetPostStepPoint();
  const CLHEP::Hep3Vector& postPos = postPoint->GetPosition();
  //Go below only in HGCal case
  if (isHGCal or isHFNose) {
    //A step never spans across boundaries: geometry or physics define the end points
    //If the step is limited by a boundary, the post-step point stands on the
    //boundary and it logically belongs to the next volume.
    if (postPoint->GetStepStatus() == fGeomBoundary && fabs(postPoint->GetMomentum().eta()) > outerHGCalEta &&
        fabs(postPoint->GetMomentum().eta()) < innerHGCalEta) {
      //Post point position is the low z edge of the new volume, or the upper for the prepoint volume.
      //So, premat - postz - posteta - postR - premattotalenergylossEtable - premattotalenergylossEfull
      //Observe the two zeros at the end which in the past where set to emCalculator.GetDEDX and
      //emCalculator.ComputeTotalDEDX but decided not to be used. Will save the structure though for the script.
      outVolumeZpositionTxt << prePoint->GetMaterial()->GetName() << " " << postPos.z() << " "
                            << postPoint->GetMomentum().eta() << " "
                            << sqrt(postPos.x() * postPos.x() + postPos.y() * postPos.y()) << " " << 0 << " " << 0
                            << std::endl;
    }

    if (postPoint->GetStepStatus() == fGeomBoundary && fabs(postPoint->GetMomentum().eta()) > outerHFnoseEta &&
        fabs(postPoint->GetMomentum().eta()) < innerHFnoseEta) {
      outVolumeZpositionTxt << prePoint->GetMaterial()->GetName() << " " << postPos.z() << " "
                            << postPoint->GetMomentum().eta() << " "
                            << sqrt(postPos.x() * postPos.x() + postPos.y() * postPos.y()) << " " << 0 << " " << 0
                            << std::endl;
    }

  }  //end of isHGCal  or HFnose if

  // check for a sensitive detector
  bool enter_sensitive = false;
  bool leave_sensitive = false;
  double cosThetaPre = 0.0;
  double cosThetaPost = 0.0;
  int level = 0;
  const G4VPhysicalVolume* sensitive = nullptr;
  GlobalPoint position;
  std::tie(sensitive, level) = GetSensitiveVolume(touchable);
  if (sensitive) {
    const G4VSolid& solid = *touchable->GetSolid(level);
    const G4AffineTransform& transform = GetTransform(touchable, level);
    G4ThreeVector pos = transform.Inverse().TransformPoint(G4ThreeVector(0., 0., 0.));
    position = GlobalPoint(pos.x() / cm, pos.y() / cm, pos.z() / cm);  // mm -> cm

    G4ThreeVector localPosPre = transform.TransformPoint(globalPosPre);
    EInside statusPre = solid.Inside(localPosPre);
    if (statusPre == kSurface) {
      enter_sensitive = true;
      G4ThreeVector globalDirPre = step->GetPreStepPoint()->GetMomentumDirection();
      G4ThreeVector localDirPre = transform.TransformAxis(globalDirPre);
      G4ThreeVector normalPre = solid.SurfaceNormal(localPosPre);
      cosThetaPre = normalPre.cosTheta(-localDirPre);
    }

    G4ThreeVector localPosPost = transform.TransformPoint(globalPosPost);
    EInside statusPost = solid.Inside(localPosPost);
    if (statusPost == kSurface) {
      leave_sensitive = true;
      G4ThreeVector globalDirPost = step->GetPostStepPoint()->GetMomentumDirection();
      G4ThreeVector localDirPost = transform.TransformAxis(globalDirPost);
      G4ThreeVector normalPost = solid.SurfaceNormal(localPosPost);
      cosThetaPost = normalPost.cosTheta(localDirPost);
    }
  }

  // update track accounting
  if (enter_sensitive) {
    if (m_track_volume != nullptr) {
      edm::LogWarning("TrackingMaterialProducer") << "Entering volume " << sensitive << " while inside volume "
                                                  << m_track_volume << ". Something is inconsistent";
      m_track.reset();
    }
    m_track_volume = sensitive;
    m_track.enterDetector(position, cosThetaPre);
  }
  m_track.step(MaterialAccountingStep(length, radiationLengths, energyLoss, globalPositionIn, globalPositionOut));
  if (leave_sensitive) {
    if (m_track_volume != sensitive) {
      edm::LogWarning("TrackingMaterialProducer") << "Leaving volume " << sensitive << " while inside volume "
                                                  << m_track_volume << ". Something is inconsistent";
      m_track.reset();
    } else
      m_track.leaveDetector(cosThetaPost);
    m_track_volume = nullptr;
  }
  if (sensitive)
    LogInfo("TrackingMaterialProducer") << "Track was near sensitive     volume " << sensitive->GetName() << std::endl;
  else
    LogInfo("TrackingMaterialProducer") << "Track was near non-sensitive volume " << touchable->GetVolume()->GetName()
                                        << std::endl;
  LogInfo("TrackingMaterialProducer") << "Step length:             " << length << " cm\n"
                                      << "globalPreStep(r,z): (" << globalPositionIn.perp() << ", "
                                      << globalPositionIn.z() << ") cm\n"
                                      << "globalPostStep(r,z): (" << globalPositionOut.perp() << ", "
                                      << globalPositionOut.z() << ") cm\n"
                                      << "position(r,z): (" << position.perp() << ", " << position.z() << ") cm\n"
                                      << "Radiation lengths:       " << radiationLengths << " \t\t(X0: " << X0
                                      << " cm)\n"
                                      << "Energy loss:             " << energyLoss << " MeV  \t(Xi: " << Xi
                                      << " MeV/cm)\n"
                                      << "Track was " << (enter_sensitive ? "entering " : "in none ")
                                      << "sensitive volume\n"
                                      << "Track was " << (leave_sensitive ? "leaving  " : "in none ")
                                      << "sensitive volume\n";
}

void TrackingMaterialProducer::update ( const EndOfTrack *  )

overrideprivatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfTrack * >.

Definition at line 331 of file TrackingMaterialProducer.cc.

References HLT_FULL_cff::track.

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                             {
  const G4Track* track = (*event)();
  if (m_primaryTracks and track->GetParentID() != 0)
    return;

  radLen_vs_eta_->Fill(track->GetMomentum().eta(), m_track.summary().radiationLengths());
  m_tracks->push_back(m_track);

  // LogInfo
  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: this track took " << m_track.steps().size()
                                      << " steps, and passed through " << m_track.detectors().size()
                                      << " sensitive detectors" << std::endl;
  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: track length:       " << m_track.summary().length()
                                      << " cm" << std::endl;
  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: radiation lengths: "
                                      << m_track.summary().radiationLengths() << std::endl;
  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: energy loss:        "
                                      << m_track.summary().energyLoss() << " MeV" << std::endl;
}

void TrackingMaterialProducer::update ( const EndOfJob *  )

overrideprivatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfJob * >.

Definition at line 114 of file TrackingMaterialProducer.cc.

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                           {
  radLen_vs_eta_->Write();
  output_file_->Close();
}

Member Data Documentation

constexpr float TrackingMaterialProducer::innerHFnoseEta = 4.

staticprivate

Definition at line 71 of file TrackingMaterialProducer.h.

constexpr float TrackingMaterialProducer::innerHGCalEta = 2.4

staticprivate

Definition at line 69 of file TrackingMaterialProducer.h.

bool TrackingMaterialProducer::isHFNose

private

Definition at line 68 of file TrackingMaterialProducer.h.

bool TrackingMaterialProducer::isHGCal

private

Definition at line 67 of file TrackingMaterialProducer.h.

double TrackingMaterialProducer::m_hgcalzfront

private

Definition at line 61 of file TrackingMaterialProducer.h.

bool TrackingMaterialProducer::m_primaryTracks

private

Definition at line 57 of file TrackingMaterialProducer.h.

std::vector<std::string> TrackingMaterialProducer::m_selectedNames

private

Definition at line 58 of file TrackingMaterialProducer.h.

std::vector<const G4LogicalVolume*> TrackingMaterialProducer::m_selectedVolumes

private

Definition at line 59 of file TrackingMaterialProducer.h.

MaterialAccountingTrack TrackingMaterialProducer::m_track

private

Definition at line 62 of file TrackingMaterialProducer.h.

const G4VPhysicalVolume* TrackingMaterialProducer::m_track_volume

private

Definition at line 63 of file TrackingMaterialProducer.h.

std::vector<MaterialAccountingTrack>* TrackingMaterialProducer::m_tracks

private

Definition at line 64 of file TrackingMaterialProducer.h.

std::string TrackingMaterialProducer::m_txtOutFile

private

Definition at line 60 of file TrackingMaterialProducer.h.

constexpr float TrackingMaterialProducer::outerHFnoseEta = 3.3

staticprivate

Definition at line 72 of file TrackingMaterialProducer.h.

constexpr float TrackingMaterialProducer::outerHGCalEta = 2.0

staticprivate

Definition at line 70 of file TrackingMaterialProducer.h.

TFile* TrackingMaterialProducer::output_file_

private

Definition at line 65 of file TrackingMaterialProducer.h.

std::ofstream TrackingMaterialProducer::outVolumeZpositionTxt

private

Definition at line 73 of file TrackingMaterialProducer.h.

TProfile* TrackingMaterialProducer::radLen_vs_eta_

private

Definition at line 66 of file TrackingMaterialProducer.h.