#include <DreamSD.h>

Inheritance diagram for DreamSD:

Public Member Functions
	DreamSD (const std::string &, const DDCompactView , const cms::DDCompactView , 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 SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, float timeSlice=1., bool ignoreTkID=false, const std::string &newcolname="")

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	newCollection (const std::string &name, edm::ParameterSet const &p)

void	PrintAll () override

G4bool	ProcessHits (G4Step step, G4TouchableHistory ) override

bool	ProcessHits (G4GFlashSpot aSpot, G4TouchableHistory ) override

void	reset () override

	~CaloSD () override

Public Member Functions inherited from SensitiveCaloDetector
	SensitiveCaloDetector (const std::string &iname, const SensitiveDetectorCatalog &clg, const std::string &newcollname="")

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 SensitiveDetectorCatalog &, bool calo, const std::string &newcollname="")

	~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 (int k=0)

CaloG4Hit *	createNewHit (const G4Step , const G4Track , int k)

virtual void	endEvent ()

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

virtual bool	filterHit (CaloG4Hit *, double)

unsigned int	findBoundaryCrossingParent (const G4Track *track, bool markParentAsSaveable=true)

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 (int k=0)

double	getResponseWt (const G4Track *, int k=0)

virtual int	getTrackID (const G4Track *)

bool	hitExists (const G4Step *, int k)

void	ignoreRejection ()

virtual void	initEvent (const BeginOfEvent *)

void	printDetectorLevels (const G4VTouchable *) const

void	processHit (const G4Step *step)

virtual void	processSecondHit (const G4Step , const G4Track )

void	resetForNewPrimary (const G4Step *)

void	setNumberCheckedHits (int val, int k=0)

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)

std::string	shortreprID (const CaloHitID &ID)

std::string	shortreprID (const CaloG4Hit *hit)

void	update (const BeginOfRun *) 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 *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	update (const ::EndOfEvent *) override

void	updateHit (CaloG4Hit *, int k)

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 &)

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

Private Attributes
double	birk1_

double	birk2_

double	birk3_

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

const cms::DDCompactView *	cpvDD4hep_

const DDCompactView *	cpvDDD_

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
std::string	collName_ [2]

CaloG4Hit *	currentHit [2]

CaloHitID	currentID [2]

std::string	detName_

float	edepositEM

float	edepositHAD

double	eminHit

double	energyCut

G4ThreeVector	entranceLocal

G4ThreeVector	entrancePoint

bool	forceSave

std::vector< std::string >	hcn_

float	incidentEnergy

double	kmaxIon

double	kmaxNeutron

double	kmaxProton

int	nHC_

G4ThreeVector	posGlobal

CaloHitID	previousID [2]

bool	suppressHeavy

double	tmaxHit

std::vector< int >	useResMap_

Detailed Description

Definition at line 20 of file DreamSD.h.

Member Typedef Documentation

◆ DimensionMap

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

private

Definition at line 38 of file DreamSD.h.

◆ Doubles

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

private

Definition at line 37 of file DreamSD.h.

Constructor & Destructor Documentation

◆ DreamSD()

DreamSD::DreamSD	(	const std::string &	name,
		const DDCompactView *	cpvDDD,
		const cms::DDCompactView *	cpvDD4hep,
		const SensitiveDetectorCatalog &	clg,
		edm::ParameterSet const &	p,
		const SimTrackManager *	manager
	)

Definition at line 32 of file DreamSD.cc.

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

     : CaloSD(name, clg, p, manager), cpvDDD_(cpvDDD), cpvDD4hep_(cpvDD4hep) {
   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);
 }

◆ ~DreamSD()

DreamSD::~DreamSD ( )

inlineoverride

Definition at line 28 of file DreamSD.h.

28 {}

Member Function Documentation

◆ cherenkovDeposit_()

double DreamSD::cherenkovDeposit_ ( const G4Step * aStep )

private

Returns the total energy due to Cherenkov radiation.

Definition at line 243 of file DreamSD.cc.

References HLT_2024v14_cff::beta, ALCARECOTkAlJpsiMuMu_cff::charge, funct::cos(), l1tPhase1JetProducer_cfi::cosPhi, Calorimetry_cff::dp, getAverageNumberOfPhotons_(), getPhotonEnergyDeposit_(), materialPropertiesTable_, phi, ElectronMcFakeValidator_cfi::Pmax, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, funct::sin(), l1tPhase1JetProducer_cfi::sinPhi, and mathSSE::sqrt().

Referenced by getEnergyDeposit().

                                                      {
   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;
 }

◆ crystalLength()

double DreamSD::crystalLength ( G4LogicalVolume * lv ) const

private

Definition at line 223 of file DreamSD.cc.

References xtalLMap_.

Referenced by curve_LY(), and getPhotonEnergyDeposit_().

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

◆ crystalWidth()

double DreamSD::crystalWidth ( G4LogicalVolume * lv ) const

private

Definition at line 232 of file DreamSD.cc.

References ApeEstimator_cff::width, and xtalLMap_.

Referenced by getPhotonEnergyDeposit_().

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

◆ curve_LY()

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

private

Definition at line 196 of file DreamSD.cc.

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

Referenced by getEnergyDeposit().

                                                       {
   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;
 }

◆ fillMap()

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

private

Definition at line 176 of file DreamSD.cc.

References mergeVDriftHistosByStation::name, DD4hep2DDDName::noNameSpace(), AlCaHLTBitMon_QueryRunRegistry::string, ApeEstimator_cff::width, and xtalLMap_.

Referenced by initMap().

                                                                         {
   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();
     std::string namex = static_cast<std::string>((*lvcite)->GetName());
     if (name == DD4hep2DDDName::noNameSpace(static_cast<std::string>(namex))) {
       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)));
 }

◆ 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 359 of file DreamSD.cc.

References HLT_2024v14_cff::beta, chAngleIntegrals_, ALCARECOTkAlJpsiMuMu_cff::charge, Calorimetry_cff::dp, metDiagnosticParameterSet_cfi::nMax, metDiagnosticParameterSet_cfi::nMin, and ElectronMcFakeValidator_cfi::Pmax.

Referenced by cherenkovDeposit_().

                                                                                    {
   const double rFact = 369.81 / (CLHEP::eV * CLHEP::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()->GetMaxEnergy();
 
   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 / CLHEP::eplus * charge / CLHEP::eplus * (dp - ge * BetaInverse * BetaInverse);
 
   edm::LogVerbatim("EcalSim") << "@SUB=getAverageNumberOfPhotons\nCAImin = " << CAImin << "\nCAImax = " << CAImax
                               << "\ndp = " << dp << ", ge = " << ge << "\nnumPhotons = " << numPhotons;
   return numPhotons;
 }

◆ getEnergyDeposit()

double DreamSD::getEnergyDeposit ( const G4Step * aStep )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 68 of file DreamSD.cc.

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

                                                     {
   // 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;
 }

◆ getPhotonEnergyDeposit_()

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

private

Returns energy deposit for a given photon.

Definition at line 505 of file DreamSD.cc.

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

Referenced by cherenkovDeposit_().

                                                                                                            {
   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;
 }

◆ initMap()

void DreamSD::initMap ( const std::string & sd )

private

Definition at line 121 of file DreamSD.cc.

References cpvDD4hep_, cpvDDD_, dd4hep_, cms::DDSolidShapeMap, fillMap(), ALCARECOTkAlBeamHalo_cff::filter, DDFilteredView::firstChild(), cms::DDFilteredView::firstChild(), mps_fire::i, k_ScaleFromDD4hepToG4, k_ScaleFromDDDToG4, DDFilteredView::logicalPart(), cms::dd::name(), mergeVDriftHistosByStation::name, cms::DDFilteredView::name(), DDFilteredView::next(), DD4hep2DDDName::noNameSpace(), cms::DDFilteredView::parameters(), cms::DDFilteredView::shape(), mkfit::Const::sol, DDLogicalPart::solid(), jetUpdater_cfi::sort, AlCaHLTBitMon_QueryRunRegistry::string, ApeEstimator_cff::width, and xtalLMap_.

Referenced by DreamSD().

                                          {
   if (dd4hep_) {
     const cms::DDFilter filter("ReadOutName", sd);
     cms::DDFilteredView fv((*cpvDD4hep_), filter);
     while (fv.firstChild()) {
       std::string name = DD4hep2DDDName::noNameSpace(static_cast<std::string>(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 {
     DDSpecificsMatchesValueFilter filter{DDValue("ReadOutName", sd, 0)};
     DDFilteredView fv((*cpvDDD_), filter);
     fv.firstChild();
     bool dodet = true;
     while (dodet) {
       const DDSolid &sol = fv.logicalPart().solid();
       std::vector<double> paras(sol.parameters());
       std::string name = static_cast<std::string>(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();
       fillMap(name, length, width);
       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
   }
 }

◆ initRun()

void DreamSD::initRun ( )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 88 of file DreamSD.cc.

References materialPropertiesTable_, setPbWO2MaterialProperties_(), and xtalLMap_.

                       {
   // 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();
   }
 }

◆ setDetUnitId()

uint32_t DreamSD::setDetUnitId ( const G4Step * aStep )

overridevirtual

Implements CaloSD.

Definition at line 107 of file DreamSD.cc.

References EcalPhiSymFlatTableProducers_cfi::id, readBothSide_, and side_.

                                                   {
   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;
 }

◆ setPbWO2MaterialProperties_()

bool DreamSD::setPbWO2MaterialProperties_ ( G4Material * aMaterial )

private

Sets material properties at run-time...

Definition at line 422 of file DreamSD.cc.

References chAngleIntegrals_, mergeVDriftHistosByStation::name, DD4hep2DDDName::noNameSpace(), AlCaHLTBitMon_QueryRunRegistry::string, and TableParser::table.

Referenced by initRun().

                                                                {
   std::string pbWO2Name("E_PbWO4");
   std::string name = static_cast<std::string>(aMaterial->GetName());
   if (DD4hep2DDDName::noNameSpace(name) != pbWO2Name) {  // 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;
   using CLHEP::eV;
   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<G4PhysicsFreeVector>(nEntries);
 
   int index = 0;
   double currentRI = RefractiveIndex[index];
   double currentPM = PhotonEnergy[index];
   double currentCAI = 0.0;
   chAngleIntegrals_.get()->PutValue(0, 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()->PutValue(index, currentPM, currentCAI);
 
     prevPM = currentPM;
     prevCAI = currentCAI;
     prevRI = currentRI;
   }
 
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("EcalSim") << "Material properties set for " << aMaterial->GetName();
 #endif
   return true;
 }

Member Data Documentation

◆ birk1_

double DreamSD::birk1_

private

Definition at line 65 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ birk2_

double DreamSD::birk2_

private

Definition at line 65 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ birk3_

double DreamSD::birk3_

private

Definition at line 65 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ chAngleIntegrals_

std::unique_ptr<G4PhysicsFreeVector> DreamSD::chAngleIntegrals_

private

Table of Cherenkov angle integrals vs photon momentum.

Definition at line 72 of file DreamSD.h.

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

◆ cpvDD4hep_

const cms::DDCompactView* DreamSD::cpvDD4hep_

private

Definition at line 62 of file DreamSD.h.

Referenced by initMap().

◆ cpvDDD_

const DDCompactView* DreamSD::cpvDDD_

private

Definition at line 61 of file DreamSD.h.

Referenced by initMap().

◆ dd4hep_

bool DreamSD::dd4hep_

private

Definition at line 64 of file DreamSD.h.

Referenced by DreamSD(), and initMap().

◆ doCherenkov_

bool DreamSD::doCherenkov_

private

Definition at line 64 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ k_ScaleFromDD4hepToG4

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

staticprivate

Definition at line 59 of file DreamSD.h.

Referenced by initMap().

◆ k_ScaleFromDDDToG4

constexpr double DreamSD::k_ScaleFromDDDToG4 = 1.0

staticprivate

Definition at line 58 of file DreamSD.h.

Referenced by initMap().

◆ materialPropertiesTable_

G4MaterialPropertiesTable* DreamSD::materialPropertiesTable_

private

Definition at line 73 of file DreamSD.h.

Referenced by cherenkovDeposit_(), and initRun().

◆ nphotons_

int DreamSD::nphotons_

private

Definition at line 75 of file DreamSD.h.

◆ readBothSide_

bool DreamSD::readBothSide_

private

Definition at line 64 of file DreamSD.h.

Referenced by DreamSD(), and setDetUnitId().

◆ side_

int DreamSD::side_

private

Definition at line 69 of file DreamSD.h.

Referenced by getEnergyDeposit(), and setDetUnitId().

◆ slopeLY_

double DreamSD::slopeLY_

private

Definition at line 66 of file DreamSD.h.

Referenced by curve_LY(), and DreamSD().

◆ useBirk_

bool DreamSD::useBirk_

private

Definition at line 64 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

◆ xtalLMap_

DimensionMap DreamSD::xtalLMap_

private

Definition at line 67 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_

◆ cpvDD4hep_

◆ cpvDDD_

◆ dd4hep_

◆ doCherenkov_

◆ k_ScaleFromDD4hepToG4

◆ k_ScaleFromDDDToG4

◆ materialPropertiesTable_

◆ nphotons_

◆ readBothSide_

◆ side_

◆ slopeLY_

◆ useBirk_

◆ xtalLMap_