#include <DreamSD.h>

Inheritance diagram for DreamSD:

Public Member Functions
	DreamSD (const std::string &, const edm::EventSetup &, const SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)

uint32_t	setDetUnitId (const G4Step *) override

	~DreamSD () override

Public Member Functions inherited from CaloSD
	CaloSD (const std::string &aSDname, const edm::EventSetup &es, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, float timeSlice=1., bool ignoreTkID=false)

void	clear () override

void	clearHits () override

void	DrawAll () override

void	EndOfEvent (G4HCofThisEvent *eventHC) override

void	fillHits (edm::PCaloHitContainer &, const std::string &) override

void	Initialize (G4HCofThisEvent *HCE) override

bool	isItFineCalo (const G4VTouchable *touch)

void	PrintAll () override

bool	ProcessHits (G4GFlashSpot aSpot, G4TouchableHistory ) override

G4bool	ProcessHits (G4Step step, G4TouchableHistory ) override

void	reset () override

	~CaloSD () override

Public Member Functions inherited from SensitiveCaloDetector
	SensitiveCaloDetector (const std::string &iname, const edm::EventSetup &es, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p)

Public Member Functions inherited from SensitiveDetector
void	EndOfEvent (G4HCofThisEvent *eventHC) override

const std::vector< std::string > &	getNames () const

void	Initialize (G4HCofThisEvent *eventHC) override

bool	isCaloSD () const

	SensitiveDetector (const std::string &iname, const edm::EventSetup &es, const SensitiveDetectorCatalog &, edm::ParameterSet const &p, bool calo)

	~SensitiveDetector () override

Public Member Functions inherited from Observer< const BeginOfRun * >
	Observer ()

void	slotForUpdate (const BeginOfRun * 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 EndOfTrack * >
	Observer ()

void	slotForUpdate (const EndOfTrack * iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const EndOfEvent * >
	Observer ()

void	slotForUpdate (const EndOfEvent * iT)

virtual	~Observer ()

Protected Member Functions
double	getEnergyDeposit (const G4Step *) override

void	initRun () override

Protected Member Functions inherited from CaloSD
bool	checkHit ()

CaloG4Hit *	createNewHit (const G4Step , const G4Track )

virtual void	endEvent ()

virtual double	EnergyCorrected (const G4Step &step, const G4Track *)

virtual bool	filterHit (CaloG4Hit *, double)

double	getAttenuation (const G4Step *aStep, double birk1, double birk2, double birk3) const

virtual uint16_t	getDepth (const G4Step *)

virtual bool	getFromLibrary (const G4Step *step)

int	getNumberOfHits ()

double	getResponseWt (const G4Track *)

virtual int	getTrackID (const G4Track *)

void	hitBookkeepingFineCalo (const G4Step step, const G4Track currentTrack, CaloG4Hit *hit)

bool	hitExists (const G4Step *)

void	ignoreRejection ()

virtual void	initEvent (const BeginOfEvent *)

void	printDetectorLevels (const G4VTouchable *) const

void	processHit (const G4Step *step)

void	resetForNewPrimary (const G4Step *)

void	setNumberCheckedHits (int val)

void	setParameterized (bool val)

G4ThreeVector	setToGlobal (const G4ThreeVector &, const G4VTouchable *) const

G4ThreeVector	setToLocal (const G4ThreeVector &, const G4VTouchable *) const

virtual int	setTrackID (const G4Step *)

void	setUseMap (bool val)

void	update (const ::EndOfEvent *) override

void	update (const BeginOfEvent *) override
	This routine will be called when the appropriate signal arrives. More...

void	update (const BeginOfRun *) override
	This routine will be called when the appropriate signal arrives. More...

void	update (const BeginOfTrack *trk) override
	This routine will be called when the appropriate signal arrives. More...

void	update (const EndOfTrack *trk) override
	This routine will be called when the appropriate signal arrives. More...

void	updateHit (CaloG4Hit *)

Protected Member Functions inherited from SensitiveDetector
TrackInformation *	cmsTrackInformation (const G4Track *aTrack)

Local3DPoint	ConvertToLocal3DPoint (const G4ThreeVector &point) const

Local3DPoint	FinalStepPosition (const G4Step *step, coordinates) const

Local3DPoint	InitialStepPosition (const G4Step *step, coordinates) const

Local3DPoint	LocalPostStepPosition (const G4Step *step) const

Local3DPoint	LocalPreStepPosition (const G4Step *step) const

void	NaNTrap (const G4Step *step) const

void	setNames (const std::vector< std::string > &)

Protected Member Functions inherited from Observer< const EndOfEvent * >
virtual void	update (const EndOfEvent *)=0
	This routine will be called when the appropriate signal arrives. More...

Private Types
typedef std::map< G4LogicalVolume *, Doubles >	DimensionMap

typedef std::pair< double, double >	Doubles

Private Member Functions
double	cherenkovDeposit_ (const G4Step *aStep)
	Returns the total energy due to Cherenkov radiation. More...

double	crystalLength (G4LogicalVolume *) const

double	crystalWidth (G4LogicalVolume *) const

double	curve_LY (const G4Step *, int)

void	fillMap (const std::string &, double, double)

double	getAverageNumberOfPhotons_ (const double charge, const double beta, const G4Material aMaterial, const G4MaterialPropertyVector rIndex)
	Returns average number of photons created by track. More...

double	getPhotonEnergyDeposit_ (const G4ParticleMomentum &p, const G4ThreeVector &x, const G4Step *aStep)
	Returns energy deposit for a given photon. More...

void	initMap (const std::string &, const edm::EventSetup &)

bool	setPbWO2MaterialProperties_ (G4Material *aMaterial)
	Sets material properties at run-time... More...

Private Attributes
double	birk1_

double	birk2_

double	birk3_

std::unique_ptr< G4PhysicsOrderedFreeVector >	chAngleIntegrals_
	Table of Cherenkov angle integrals vs photon momentum. More...

bool	dd4hep_

bool	doCherenkov_

G4MaterialPropertiesTable *	materialPropertiesTable_

int	nphotons_

bool	readBothSide_

int	side_

double	slopeLY_

bool	useBirk_

DimensionMap	xtalLMap_

Static Private Attributes
static constexpr double	k_ScaleFromDD4HepToG4 = 1.0 / dd4hep::mm

static constexpr double	k_ScaleFromDDDToG4 = 1.0

Additional Inherited Members
Protected Types inherited from SensitiveDetector
enum	coordinates { WorldCoordinates, LocalCoordinates }

Static Protected Member Functions inherited from CaloSD
static std::string	printableDecayChain (const std::vector< unsigned int > &decayChain)

Protected Attributes inherited from CaloSD
CaloG4Hit *	currentHit

CaloHitID	currentID

float	edepositEM

float	edepositHAD

double	eminHit

double	energyCut

G4ThreeVector	entranceLocal

G4ThreeVector	entrancePoint

bool	forceSave

float	incidentEnergy

double	kmaxIon

double	kmaxNeutron

double	kmaxProton

G4ThreeVector	posGlobal

CaloHitID	previousID

bool	suppressHeavy

double	tmaxHit

Detailed Description

Definition at line 17 of file DreamSD.h.

Member Typedef Documentation

◆ DimensionMap

typedef std::map<G4LogicalVolume *, Doubles> DreamSD::DimensionMap

private

Definition at line 34 of file DreamSD.h.

◆ Doubles

typedef std::pair<double, double> DreamSD::Doubles

private

Definition at line 33 of file DreamSD.h.

Constructor & Destructor Documentation

◆ DreamSD()

DreamSD::DreamSD	(	const std::string &	name,
		const edm::EventSetup &	es,
		const SensitiveDetectorCatalog &	clg,
		edm::ParameterSet const &	p,
		const SimTrackManager *	manager
	)

Definition at line 37 of file DreamSD.cc.

     : CaloSD(name, es, clg, p, manager) {
   edm::ParameterSet m_EC = p.getParameter<edm::ParameterSet>("ECalSD");
   useBirk_ = m_EC.getParameter<bool>("UseBirkLaw");
   doCherenkov_ = m_EC.getParameter<bool>("doCherenkov");
   birk1_ = m_EC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
   birk2_ = m_EC.getParameter<double>("BirkC2");
   birk3_ = m_EC.getParameter<double>("BirkC3");
   slopeLY_ = m_EC.getParameter<double>("SlopeLightYield");
   readBothSide_ = m_EC.getUntrackedParameter<bool>("ReadBothSide", false);
   dd4hep_ = p.getParameter<bool>("g4GeometryDD4hepSource");
  
   chAngleIntegrals_.reset(nullptr);
  
   edm::LogVerbatim("EcalSim") << "Constructing a DreamSD  with name " << GetName()
                               << "\nDreamSD:: Use of Birks law is set to      " << useBirk_
                               << "  with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_
                               << "\n          Slope for Light yield is set to " << slopeLY_
                               << "\n          Parameterization of Cherenkov is set to " << doCherenkov_
                               << ", readout both sides is " << readBothSide_ << " and dd4hep flag " << dd4hep_;
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << GetName() << " initialized";
   const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
   unsigned int k(0);
   for (auto lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite, ++k)
     edm::LogVerbatim("EcalSim") << "Volume[" << k << "] " << (*lvcite)->GetName();
 #endif
   initMap(name, es);
 }

References birk1_, birk2_, birk3_, chAngleIntegrals_, dd4hep_, doCherenkov_, g, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), initMap(), dqmdumpme::k, MeV, Skims_PA_cff::name, AlCaHLTBitMon_ParallelJobs::p, readBothSide_, slopeLY_, and useBirk_.

◆ ~DreamSD()

DreamSD::~DreamSD ( )

inlineoverride

Definition at line 24 of file DreamSD.h.

24 {}

Member Function Documentation

◆ cherenkovDeposit_()

double DreamSD::cherenkovDeposit_ ( const G4Step * aStep )

private

Returns the total energy due to Cherenkov radiation.

Definition at line 260 of file DreamSD.cc.

                                                      {
   double cherenkovEnergy = 0;
   if (!materialPropertiesTable_)
     return cherenkovEnergy;
   G4Material *material = aStep->GetTrack()->GetMaterial();
  
   // Retrieve refractive index
   G4MaterialPropertyVector *Rindex = materialPropertiesTable_->GetProperty("RINDEX");
   if (Rindex == nullptr) {
     edm::LogWarning("EcalSim") << "Couldn't retrieve refractive index";
     return cherenkovEnergy;
   }
  
   // V.Ivanchenko - temporary close log output for 9.5
   // Material refraction properties
   int Rlength = Rindex->GetVectorLength() - 1;
   double Pmin = Rindex->Energy(0);
   double Pmax = Rindex->Energy(Rlength);
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "Material properties: \n  Pmin = " << Pmin << "  Pmax = " << Pmax;
 #endif
   // Get particle properties
   const G4StepPoint *pPreStepPoint = aStep->GetPreStepPoint();
   const G4StepPoint *pPostStepPoint = aStep->GetPostStepPoint();
   const G4ThreeVector &x0 = pPreStepPoint->GetPosition();
   G4ThreeVector p0 = aStep->GetDeltaPosition().unit();
   const G4DynamicParticle *aParticle = aStep->GetTrack()->GetDynamicParticle();
   const double charge = aParticle->GetDefinition()->GetPDGCharge();
   // beta is averaged over step
   double beta = 0.5 * (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta());
   double BetaInverse = 1.0 / beta;
  
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "Particle properties: \n  charge = " << charge << "  beta   = " << beta;
 #endif
  
   // Now get number of photons generated in this step
   double meanNumberOfPhotons = getAverageNumberOfPhotons_(charge, beta, material, Rindex);
   if (meanNumberOfPhotons <= 0.0) {  // Don't do anything
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("EcalSim") << "Mean number of photons is zero: " << meanNumberOfPhotons << ", stopping here";
 #endif
     return cherenkovEnergy;
   }
  
   // number of photons is in unit of Geant4...
   meanNumberOfPhotons *= aStep->GetStepLength();
  
   // Now get a poisson distribution
   int numPhotons = static_cast<int>(G4Poisson(meanNumberOfPhotons));
   // edm::LogVerbatim("EcalSim") << "Number of photons = " << numPhotons;
   if (numPhotons <= 0) {
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("EcalSim") << "Poission number of photons is zero: " << numPhotons << ", stopping here";
 #endif
     return cherenkovEnergy;
   }
  
   // Material refraction properties
   double dp = Pmax - Pmin;
   double maxCos = BetaInverse / (*Rindex)[Rlength];
   double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
  
   // Finally: get contribution of each photon
   for (int iPhoton = 0; iPhoton < numPhotons; ++iPhoton) {
     // Determine photon momentum
     double randomNumber;
     double sampledMomentum, sampledRI;
     double cosTheta, sin2Theta;
  
     // sample a momentum (not sure why this is needed!)
     do {
       randomNumber = G4UniformRand();
       sampledMomentum = Pmin + randomNumber * dp;
       sampledRI = Rindex->Value(sampledMomentum);
       cosTheta = BetaInverse / sampledRI;
  
       sin2Theta = (1.0 - cosTheta) * (1.0 + cosTheta);
       randomNumber = G4UniformRand();
  
     } while (randomNumber * maxSin2 > sin2Theta);
  
     // Generate random position of photon on cone surface
     // defined by Theta
     randomNumber = G4UniformRand();
  
     double phi = twopi * randomNumber;
     double sinPhi = sin(phi);
     double cosPhi = cos(phi);
  
     // Create photon momentum direction vector
     // The momentum direction is still w.r.t. the coordinate system where the
     // primary particle direction is aligned with the z axis
     double sinTheta = sqrt(sin2Theta);
     double px = sinTheta * cosPhi;
     double py = sinTheta * sinPhi;
     double pz = cosTheta;
     G4ThreeVector photonDirection(px, py, pz);
  
     // Rotate momentum direction back to global (crystal) reference system
     photonDirection.rotateUz(p0);
  
     // Create photon position and momentum
     randomNumber = G4UniformRand();
     G4ThreeVector photonPosition = x0 + randomNumber * aStep->GetDeltaPosition();
     G4ThreeVector photonMomentum = sampledMomentum * photonDirection;
  
     // Collect energy on APD
     cherenkovEnergy += getPhotonEnergyDeposit_(photonMomentum, photonPosition, aStep);
   }
   return cherenkovEnergy;
 }

References HLT_FULL_cff::beta, ALCARECOTkAlJpsiMuMu_cff::charge, funct::cos(), Calorimetry_cff::dp, getAverageNumberOfPhotons_(), getPhotonEnergyDeposit_(), materialPropertiesTable_, phi, dataAnalyzerFineBiningParameters_cff::Pmax, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, funct::sin(), and mathSSE::sqrt().

Referenced by getEnergyDeposit().

◆ crystalLength()

double DreamSD::crystalLength ( G4LogicalVolume * lv ) const

private

Definition at line 240 of file DreamSD.cc.

                                                        {
   double length = -1.;
   DimensionMap::const_iterator ite = xtalLMap_.find(lv);
   if (ite != xtalLMap_.end())
     length = ite->second.first;
   return length;
 }

References xtalLMap_.

Referenced by curve_LY(), and getPhotonEnergyDeposit_().

◆ crystalWidth()

double DreamSD::crystalWidth ( G4LogicalVolume * lv ) const

private

Definition at line 249 of file DreamSD.cc.

                                                       {
   double width = -1.;
   DimensionMap::const_iterator ite = xtalLMap_.find(lv);
   if (ite != xtalLMap_.end())
     width = ite->second.second;
   return width;
 }

References ApeEstimator_cff::width, and xtalLMap_.

Referenced by getPhotonEnergyDeposit_().

◆ curve_LY()

double DreamSD::curve_LY	(	const G4Step *	aStep,
		int	flag
	)

private

Definition at line 213 of file DreamSD.cc.

                                                       {
   auto const stepPoint = aStep->GetPreStepPoint();
   auto const lv = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
   G4String nameVolume = lv->GetName();
  
   double weight = 1.;
   G4ThreeVector localPoint = setToLocal(stepPoint->GetPosition(), stepPoint->GetTouchable());
   double crlength = crystalLength(lv);
   double localz = localPoint.x();
   double dapd = 0.5 * crlength - flag * localz;  // Distance from closest APD
   if (dapd >= -0.1 || dapd <= crlength + 0.1) {
     if (dapd <= 100.)
       weight = 1.0 + slopeLY_ - dapd * 0.01 * slopeLY_;
   } else {
     edm::LogWarning("EcalSim") << "DreamSD: light coll curve : wrong distance "
                                << "to APD " << dapd << " crlength = " << crlength << " crystal name = " << nameVolume
                                << " z of localPoint = " << localz << " take weight = " << weight;
   }
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "DreamSD, light coll curve : " << dapd << " crlength = " << crlength
                               << " crystal name = " << nameVolume << " z of localPoint = " << localz
                               << " take weight = " << weight;
 #endif
   return weight;
 }

References crystalLength(), RemoveAddSevLevel::flag, CaloSD::setToLocal(), slopeLY_, and mps_merge::weight.

Referenced by getEnergyDeposit().

◆ fillMap()

void DreamSD::fillMap	(	const std::string &	name,
		double	length,
		double	width
	)

private

Definition at line 194 of file DreamSD.cc.

                                                                         {
   const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
   edm::LogVerbatim("EcalSim") << "LV Store with " << lvs->size() << " elements";
   G4LogicalVolume *lv = nullptr;
   for (auto lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
     edm::LogVerbatim("EcalSim") << name << " vs " << (*lvcite)->GetName();
     if ((*lvcite)->GetName() == static_cast<G4String>(name)) {
       lv = (*lvcite);
       break;
     }
   }
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "DreamSD::fillMap (for " << name << " Logical Volume " << lv << " Length " << length
                               << " Width " << width;
 #endif
   xtalLMap_.insert(std::pair<G4LogicalVolume *, Doubles>(lv, Doubles(length, width)));
 }

References Skims_PA_cff::name, ApeEstimator_cff::width, and xtalLMap_.

Referenced by initMap().

◆ getAverageNumberOfPhotons_()

double DreamSD::getAverageNumberOfPhotons_	(	const double	charge,
		const double	beta,
		const G4Material *	aMaterial,
		const G4MaterialPropertyVector *	rIndex
	)

private

Returns average number of photons created by track.

Definition at line 376 of file DreamSD.cc.

                                                                                    {
   const double rFact = 369.81 / (eV * cm);
  
   if (beta <= 0.0)
     return 0.0;
  
   double BetaInverse = 1. / beta;
  
   // Vectors used in computation of Cerenkov Angle Integral:
   //         - Refraction Indices for the current material
   //        - new G4PhysicsOrderedFreeVector allocated to hold CAI's
  
   // Min and Max photon momenta
   int Rlength = Rindex->GetVectorLength() - 1;
   double Pmin = Rindex->Energy(0);
   double Pmax = Rindex->Energy(Rlength);
  
   // Min and Max Refraction Indices
   double nMin = (*Rindex)[0];
   double nMax = (*Rindex)[Rlength];
  
   // Max Cerenkov Angle Integral
   double CAImax = chAngleIntegrals_.get()->GetMaxValue();
  
   double dp = 0., ge = 0., CAImin = 0.;
  
   // If n(Pmax) < 1/Beta -- no photons generated
   if (nMax < BetaInverse) {
   }
  
   // otherwise if n(Pmin) >= 1/Beta -- photons generated
   else if (nMin > BetaInverse) {
     dp = Pmax - Pmin;
     ge = CAImax;
   }
   // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
   // we need to find a P such that the value of n(P) == 1/Beta.
   // Interpolation is performed by the GetPhotonEnergy() and
   // GetProperty() methods of the G4MaterialPropertiesTable and
   // the GetValue() method of G4PhysicsVector.
   else {
     Pmin = Rindex->Value(BetaInverse);
     dp = Pmax - Pmin;
     // need boolean for current implementation of G4PhysicsVector
     // ==> being phased out
     double CAImin = chAngleIntegrals_->Value(Pmin);
     ge = CAImax - CAImin;
   }
  
   // Calculate number of photons
   double numPhotons = rFact * charge / eplus * charge / eplus * (dp - ge * BetaInverse * BetaInverse);
  
   edm::LogVerbatim("EcalSim") << "@SUB=getAverageNumberOfPhotons\nCAImin = " << CAImin << "\nCAImax = " << CAImax
                               << "\ndp = " << dp << ", ge = " << ge << "\nnumPhotons = " << numPhotons;
   return numPhotons;
 }

References HLT_FULL_cff::beta, chAngleIntegrals_, ALCARECOTkAlJpsiMuMu_cff::charge, Calorimetry_cff::dp, metDiagnosticParameterSet_cfi::nMax, metDiagnosticParameterSet_cfi::nMin, and dataAnalyzerFineBiningParameters_cff::Pmax.

Referenced by cherenkovDeposit_().

◆ getEnergyDeposit()

double DreamSD::getEnergyDeposit ( const G4Step * aStep )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 72 of file DreamSD.cc.

                                                     {
   // take into account light collection curve for crystals
   double weight = curve_LY(aStep, side_);
   if (useBirk_)
     weight *= getAttenuation(aStep, birk1_, birk2_, birk3_);
   double edep = aStep->GetTotalEnergyDeposit() * weight;
  
   // Get Cerenkov contribution
   if (doCherenkov_) {
     edep += cherenkovDeposit_(aStep);
   }
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "DreamSD:: " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() << " Side "
                               << side_ << " Light Collection Efficiency " << weight << " Weighted Energy Deposit "
                               << edep / CLHEP::MeV << " MeV";
 #endif
   return edep;
 }

References birk1_, birk2_, birk3_, cherenkovDeposit_(), curve_LY(), doCherenkov_, CaloSD::getAttenuation(), MeV, side_, useBirk_, and mps_merge::weight.

◆ getPhotonEnergyDeposit_()

double DreamSD::getPhotonEnergyDeposit_	(	const G4ParticleMomentum &	p,
		const G4ThreeVector &	x,
		const G4Step *	aStep
	)

private

Returns energy deposit for a given photon.

Definition at line 520 of file DreamSD.cc.

                                                                                                            {
   double energy = 0;
  
   // Crystal dimensions
  
   // edm::LogVerbatim("EcalSim") << p << x;
  
   // 1. Check if this photon goes straight to the APD:
   //    - assume that APD is at x=xtalLength/2.0
   //    - extrapolate from x=x0 to x=xtalLength/2.0 using momentum in x-y
  
   G4StepPoint *stepPoint = aStep->GetPreStepPoint();
   G4LogicalVolume *lv = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
   G4String nameVolume = lv->GetName();
  
   double crlength = crystalLength(lv);
   double crwidth = crystalWidth(lv);
   double dapd = 0.5 * crlength - x.x();  // Distance from closest APD
   double y = p.y() / p.x() * dapd;
  
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "Distance to APD: " << dapd << " - y at APD: " << y;
 #endif
   // Not straight: compute probability
   if (std::abs(y) > crwidth * 0.5) {
   }
  
   // 2. Retrieve efficiency for this wavelength (in nm, from MeV)
   double waveLength = p.mag() * 1.239e8;
  
   energy = p.mag() * PMTResponse::getEfficiency(waveLength);
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "Wavelength: " << waveLength << " - Energy: " << energy;
 #endif
   return energy;
 }

References funct::abs(), crystalLength(), crystalWidth(), HCALHighEnergyHPDFilter_cfi::energy, PMTResponse::getEfficiency(), AlCaHLTBitMon_ParallelJobs::p, x, and y.

Referenced by cherenkovDeposit_().

◆ initMap()

void DreamSD::initMap	(	const std::string &	sd,
		const edm::EventSetup &	es
	)

private

Definition at line 125 of file DreamSD.cc.

                                                                   {
   if (dd4hep_) {
     edm::ESTransientHandle<cms::DDCompactView> cpv;
     es.get<IdealGeometryRecord>().get(cpv);
     const cms::DDFilter filter("ReadOutName", sd);
     cms::DDFilteredView fv((*cpv), filter);
     while (fv.firstChild()) {
       std::string name = static_cast<std::string>(dd4hep::dd::noNamespace(fv.name()));
       std::vector<double> paras(fv.parameters());
 #ifdef EDM_ML_DEBUG
       edm::LogVerbatim("EcalSim") << "DreamSD::initMap (for " << sd << "): Solid " << name << " Shape "
                                   << cms::dd::name(cms::DDSolidShapeMap, fv.shape()) << " Parameter 0 = " << paras[0];
 #endif
       // Set length to be the largest size, width the smallest
       std::sort(paras.begin(), paras.end());
       double length = 2.0 * k_ScaleFromDD4HepToG4 * paras.back();
       double width = 2.0 * k_ScaleFromDD4HepToG4 * paras.front();
       fillMap(name, length, width);
     }
   } else {
     edm::ESTransientHandle<DDCompactView> cpv;
     es.get<IdealGeometryRecord>().get(cpv);
     DDSpecificsMatchesValueFilter filter{DDValue("ReadOutName", sd, 0)};
     DDFilteredView fv((*cpv), filter);
     fv.firstChild();
     bool dodet = true;
     const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
     while (dodet) {
       const DDSolid &sol = fv.logicalPart().solid();
       std::vector<double> paras(sol.parameters());
       G4String name = sol.name().name();
 #ifdef EDM_ML_DEBUG
       edm::LogVerbatim("EcalSim") << "DreamSD::initMap (for " << sd << "): Solid " << name << " Shape " << sol.shape()
                                   << " Parameter 0 = " << paras[0];
 #endif
       // Set length to be the largest size, width the smallest
       std::sort(paras.begin(), paras.end());
       double length = 2.0 * k_ScaleFromDDDToG4 * paras.back();
       double width = 2.0 * k_ScaleFromDDDToG4 * paras.front();
       G4LogicalVolume *lv = nullptr;
       for (auto lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
         if ((*lvcite)->GetName() == name) {
           lv = (*lvcite);
           break;
         }
       xtalLMap_.insert(std::pair<G4LogicalVolume *, Doubles>(lv, Doubles(length, width)));
 #ifdef EDM_ML_DEBUG
       edm::LogVerbatim("EcalSim") << "DreamSD " << name << ":" << lv << ":" << length << ":" << width;
 #endif
       dodet = fv.next();
     }
   }
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "DreamSD: Length Table for ReadOutName = " << sd << ":";
 #endif
   DimensionMap::const_iterator ite = xtalLMap_.begin();
   int i = 0;
   for (; ite != xtalLMap_.end(); ite++, i++) {
     G4String name = "Unknown";
     if (ite->first != nullptr)
       name = (ite->first)->GetName();
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("EcalSim") << " " << i << " " << ite->first << " " << name << " L = " << ite->second.first
                                 << " W = " << ite->second.second;
 #endif
   }
 }

References dd4hep_, cms::DDSolidShapeMap, fillMap(), ALCARECOTkAlBeamHalo_cff::filter, DDFilteredView::firstChild(), cms::DDFilteredView::firstChild(), edm::EventSetup::get(), get, mps_fire::i, k_ScaleFromDD4HepToG4, k_ScaleFromDDDToG4, DDFilteredView::logicalPart(), Skims_PA_cff::name, cms::dd::name(), DDName::name(), DDBase< N, C >::name(), cms::DDFilteredView::name(), DDFilteredView::next(), DDSolid::parameters(), cms::DDFilteredView::parameters(), sd, DDSolid::shape(), cms::DDFilteredView::shape(), DDLogicalPart::solid(), jetUpdater_cfi::sort, AlCaHLTBitMon_QueryRunRegistry::string, ApeEstimator_cff::width, and xtalLMap_.

Referenced by DreamSD().

◆ initRun()

void DreamSD::initRun ( )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 92 of file DreamSD.cc.

                       {
   // Get the material and set properties if needed
   DimensionMap::const_iterator ite = xtalLMap_.begin();
   const G4LogicalVolume *lv = (ite->first);
   G4Material *material = lv->GetMaterial();
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "DreamSD::initRun: Initializes for material " << material->GetName() << " in "
                               << lv->GetName();
 #endif
   materialPropertiesTable_ = material->GetMaterialPropertiesTable();
   if (!materialPropertiesTable_) {
     if (!setPbWO2MaterialProperties_(material)) {
       edm::LogWarning("EcalSim") << "Couldn't retrieve material properties table\n Material = " << material->GetName();
     }
     materialPropertiesTable_ = material->GetMaterialPropertiesTable();
   }
 }

References materialPropertiesTable_, setPbWO2MaterialProperties_(), and xtalLMap_.

◆ setDetUnitId()

uint32_t DreamSD::setDetUnitId ( const G4Step * aStep )

overridevirtual

Implements CaloSD.

Definition at line 111 of file DreamSD.cc.

                                                   {
   const G4VTouchable *touch = aStep->GetPreStepPoint()->GetTouchable();
   uint32_t id = (touch->GetReplicaNumber(1)) * 10 + (touch->GetReplicaNumber(0));
   side_ = readBothSide_ ? -1 : 1;
   if (side_ < 0) {
     ++id;
   }
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "DreamSD:: ID " << id;
 #endif
   return id;
 }

References triggerObjects_cff::id, readBothSide_, and side_.

◆ setPbWO2MaterialProperties_()

bool DreamSD::setPbWO2MaterialProperties_ ( G4Material * aMaterial )

private

Sets material properties at run-time...

Definition at line 439 of file DreamSD.cc.

                                                                {
   std::string pbWO2Name("E_PbWO4");
   if (pbWO2Name != aMaterial->GetName()) {  // Wrong material!
     edm::LogWarning("EcalSim") << "This is not the right material: "
                                << "expecting " << pbWO2Name << ", got " << aMaterial->GetName();
     return false;
   }
  
   G4MaterialPropertiesTable *table = new G4MaterialPropertiesTable();
  
   // Refractive index as a function of photon momentum
   // FIXME: Should somehow put that in the configuration
   const int nEntries = 14;
   double PhotonEnergy[nEntries] = {1.7712 * eV,
                                    1.8368 * eV,
                                    1.90745 * eV,
                                    1.98375 * eV,
                                    2.0664 * eV,
                                    2.15625 * eV,
                                    2.25426 * eV,
                                    2.3616 * eV,
                                    2.47968 * eV,
                                    2.61019 * eV,
                                    2.75521 * eV,
                                    2.91728 * eV,
                                    3.09961 * eV,
                                    3.30625 * eV};
   double RefractiveIndex[nEntries] = {2.17728,
                                       2.18025,
                                       2.18357,
                                       2.18753,
                                       2.19285,
                                       2.19813,
                                       2.20441,
                                       2.21337,
                                       2.22328,
                                       2.23619,
                                       2.25203,
                                       2.27381,
                                       2.30282,
                                       2.34666};
  
   table->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex, nEntries);
   aMaterial->SetMaterialPropertiesTable(table);  // FIXME: could this leak? What does G4 do?
  
   // Calculate Cherenkov angle integrals:
   // This is an ad-hoc solution (we hold it in the class, not in the material)
   chAngleIntegrals_ = std::make_unique<G4PhysicsOrderedFreeVector>();
  
   int index = 0;
   double currentRI = RefractiveIndex[index];
   double currentPM = PhotonEnergy[index];
   double currentCAI = 0.0;
   chAngleIntegrals_.get()->InsertValues(currentPM, currentCAI);
   double prevPM = currentPM;
   double prevCAI = currentCAI;
   double prevRI = currentRI;
   while (++index < nEntries) {
     currentRI = RefractiveIndex[index];
     currentPM = PhotonEnergy[index];
     currentCAI = 0.5 * (1.0 / (prevRI * prevRI) + 1.0 / (currentRI * currentRI));
     currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
  
     chAngleIntegrals_.get()->InsertValues(currentPM, currentCAI);
  
     prevPM = currentPM;
     prevCAI = currentCAI;
     prevRI = currentRI;
   }
  
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "Material properties set for " << aMaterial->GetName();
 #endif
   return true;
 }

References chAngleIntegrals_, AlCaHLTBitMon_QueryRunRegistry::string, and TableParser::table.

Referenced by initRun().

Member Data Documentation

◆ birk1_

double DreamSD::birk1_

private

Definition at line 58 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ birk2_

double DreamSD::birk2_

private

Definition at line 58 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ birk3_

double DreamSD::birk3_

private

Definition at line 58 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ chAngleIntegrals_

std::unique_ptr<G4PhysicsOrderedFreeVector> DreamSD::chAngleIntegrals_

private

Table of Cherenkov angle integrals vs photon momentum.

Definition at line 65 of file DreamSD.h.

Referenced by DreamSD(), getAverageNumberOfPhotons_(), and setPbWO2MaterialProperties_().

◆ dd4hep_

bool DreamSD::dd4hep_

private

Definition at line 57 of file DreamSD.h.

Referenced by DreamSD(), and initMap().

◆ doCherenkov_

bool DreamSD::doCherenkov_

private

Definition at line 57 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ k_ScaleFromDD4HepToG4

constexpr double DreamSD::k_ScaleFromDD4HepToG4 = 1.0 / dd4hep::mm

staticconstexprprivate

Definition at line 55 of file DreamSD.h.

Referenced by initMap().

◆ k_ScaleFromDDDToG4

constexpr double DreamSD::k_ScaleFromDDDToG4 = 1.0

staticconstexprprivate

Definition at line 54 of file DreamSD.h.

Referenced by initMap().

◆ materialPropertiesTable_

G4MaterialPropertiesTable* DreamSD::materialPropertiesTable_

private

Definition at line 66 of file DreamSD.h.

Referenced by cherenkovDeposit_(), and initRun().

◆ nphotons_

int DreamSD::nphotons_

private

Definition at line 68 of file DreamSD.h.

◆ readBothSide_

bool DreamSD::readBothSide_

private

Definition at line 57 of file DreamSD.h.

Referenced by DreamSD(), and setDetUnitId().

◆ side_

int DreamSD::side_

private

Definition at line 62 of file DreamSD.h.

Referenced by getEnergyDeposit(), and setDetUnitId().

◆ slopeLY_

double DreamSD::slopeLY_

private

Definition at line 59 of file DreamSD.h.

Referenced by curve_LY(), and DreamSD().

◆ useBirk_

bool DreamSD::useBirk_

private

Definition at line 57 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ xtalLMap_

DimensionMap DreamSD::xtalLMap_

private

Definition at line 60 of file DreamSD.h.

Referenced by crystalLength(), crystalWidth(), fillMap(), initMap(), and initRun().

Public Member Functions

Protected Member Functions

Private Types

Private Member Functions

Private Attributes

Static Private Attributes

Additional Inherited Members

Detailed Description

Member Typedef Documentation

◆ DimensionMap

◆ Doubles

Constructor & Destructor Documentation

◆ DreamSD()

◆ ~DreamSD()

Member Function Documentation

◆ cherenkovDeposit_()

◆ crystalLength()

◆ crystalWidth()

◆ curve_LY()

◆ fillMap()

◆ getAverageNumberOfPhotons_()

◆ getEnergyDeposit()

◆ getPhotonEnergyDeposit_()

◆ initMap()

◆ initRun()

◆ setDetUnitId()

◆ setPbWO2MaterialProperties_()

Member Data Documentation

◆ birk1_

◆ birk2_

◆ birk3_

◆ chAngleIntegrals_

◆ dd4hep_

◆ doCherenkov_

◆ k_ScaleFromDD4HepToG4

◆ k_ScaleFromDDDToG4

◆ materialPropertiesTable_

◆ nphotons_

◆ readBothSide_

◆ side_

◆ slopeLY_

◆ useBirk_

◆ xtalLMap_