#include <DreamSD.h>

Inheritance diagram for DreamSD:

Public Member Functions
	DreamSD (const std::string &, const 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 DDCompactView &cpv, 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

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 DDCompactView &cpv, 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 DDCompactView &cpv, 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 *)

virtual void	endEvent ()

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

bool	hitExists (const G4Step *)

virtual void	initEvent (const BeginOfEvent *)

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

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)

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

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	doCherenkov_

G4MaterialPropertiesTable *	materialPropertiesTable

int	nphotons_

bool	readBothSide_

int	side

double	slopeLY

bool	useBirk

DimensionMap	xtalLMap

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

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 15 of file DreamSD.h.

Member Typedef Documentation

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

private

Definition at line 32 of file DreamSD.h.

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

private

Definition at line 31 of file DreamSD.h.

Constructor & Destructor Documentation

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

Definition at line 27 of file DreamSD.cc.

References birk1, birk2, birk3, chAngleIntegrals_, doCherenkov_, g, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), initMap(), MeV, readBothSide_, slopeLY, and useBirk.

     : CaloSD(name, cpv, 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") * (g / (MeV * 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);
 
   chAngleIntegrals_.reset(nullptr);
 
   edm::LogInfo("EcalSim") << "Constructing a DreamSD  with name " << GetName() << "\n"
                           << "DreamSD:: 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_
                           << " and readout both sides is " << readBothSide_;
 
   initMap(name, cpv);
 }

DreamSD::~DreamSD ( )

inlineoverride

Definition at line 22 of file DreamSD.h.

References getEnergyDeposit(), initRun(), and setDetUnitId().

22 {}

Member Function Documentation

double DreamSD::cherenkovDeposit_ ( const G4Step * aStep )

private

Returns the total energy due to Cherenkov radiation.

Definition at line 191 of file DreamSD.cc.

References pfBoostedDoubleSVAK8TagInfos_cfi::beta, ALCARECOTkAlJpsiMuMu_cff::charge, funct::cos(), getAverageNumberOfPhotons_(), getPhotonEnergyDeposit_(), LogDebug, materialPropertiesTable, phi, dataAnalyzerFineBiningParameters_cff::Pmax, funct::sin(), 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);
   LogDebug("EcalSim") << "Material properties: "
                       << "\n"
                       << "  Pmin = " << Pmin << "  Pmax = " << Pmax;
 
   // 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;
 
   LogDebug("EcalSim") << "Particle properties: "
                       << "\n"
                       << "  charge = " << charge << "  beta   = " << beta;
 
   // Now get number of photons generated in this step
   double meanNumberOfPhotons = getAverageNumberOfPhotons_(charge, beta, material, Rindex);
   if (meanNumberOfPhotons <= 0.0) {  // Don't do anything
     LogDebug("EcalSim") << "Mean number of photons is zero: " << meanNumberOfPhotons << ", stopping here";
     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) {
     LogDebug("EcalSim") << "Poission number of photons is zero: " << numPhotons << ", stopping here";
     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;
 }

double DreamSD::crystalLength ( G4LogicalVolume * lv ) const

private

Definition at line 171 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;
 }

double DreamSD::crystalWidth ( G4LogicalVolume * lv ) const

private

Definition at line 180 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;
 }

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

private

Definition at line 146 of file DreamSD.cc.

References crystalLength(), LogDebug, 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;
   }
   LogDebug("EcalSim") << "DreamSD, light coll curve : " << dapd << " crlength = " << crlength
                       << " crystal name = " << nameVolume << " z of localPoint = " << localz
                       << " take weight = " << weight;
   return weight;
 }

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

References pfBoostedDoubleSVAK8TagInfos_cfi::beta, chAngleIntegrals_, LogDebug, and dataAnalyzerFineBiningParameters_cff::Pmax.

Referenced by cherenkovDeposit_().

                                                                                    {
   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);
 
   LogDebug("EcalSim") << "@SUB=getAverageNumberOfPhotons"
                       << "CAImin = " << CAImin << "\n"
                       << "CAImax = " << CAImax << "\n"
                       << "dp = " << dp << ", ge = " << ge << "\n"
                       << "numPhotons = " << numPhotons;
 
   return numPhotons;
 }

double DreamSD::getEnergyDeposit ( const G4Step * aStep )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 55 of file DreamSD.cc.

References birk1, birk2, birk3, cherenkovDeposit_(), curve_LY(), doCherenkov_, CaloSD::getAttenuation(), LogDebug, MeV, side, useBirk, and mps_merge::weight.

Referenced by ~DreamSD().

                                                     {
   // 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);
   }
   LogDebug("EcalSim") << "DreamSD:: " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() << " Side " << side
                       << " Light Collection Efficiency " << weight << " Weighted Energy Deposit " << edep / MeV
                       << " MeV";
 
   return edep;
 }

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

private

Returns energy deposit for a given photon.

Definition at line 451 of file DreamSD.cc.

References funct::abs(), crystalLength(), crystalWidth(), PMTResponse::getEfficiency(), LogDebug, 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;
 
   LogDebug("EcalSim") << "Distance to APD: " << dapd << " - y at APD: " << y;
 
   // 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);
 
   LogDebug("EcalSim") << "Wavelength: " << waveLength << " - Energy: " << energy;
 
   return energy;
 }

void DreamSD::initMap	(	const std::string &	sd,
		const DDCompactView &	cpv
	)

private

Definition at line 104 of file DreamSD.cc.

References ALCARECOTkAlBeamHalo_cff::filter, DDFilteredView::firstChild(), mps_fire::i, LogDebug, DDFilteredView::logicalPart(), DDName::name(), dataset::name, DDBase< N, C >::name(), DDFilteredView::next(), DDSolid::parameters(), DDSolid::shape(), DDLogicalPart::solid(), ApeEstimator_cff::width, and xtalLMap.

Referenced by DreamSD().

                                                                    {
   G4String attribute = "ReadOutName";
   DDSpecificsMatchesValueFilter filter{DDValue(attribute, sd, 0)};
   DDFilteredView fv(cpv, filter);
   fv.firstChild();
 
   const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
   std::vector<G4LogicalVolume *>::const_iterator lvcite;
   bool dodet = true;
   while (dodet) {
     const DDSolid &sol = fv.logicalPart().solid();
     std::vector<double> paras(sol.parameters());
     G4String name = sol.name().name();
     G4LogicalVolume *lv = nullptr;
     for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
       if ((*lvcite)->GetName() == name) {
         lv = (*lvcite);
         break;
       }
     LogDebug("EcalSim") << "DreamSD::initMap (for " << sd << "): Solid " << name << " Shape " << sol.shape()
                         << " Parameter 0 = " << paras[0] << " Logical Volume " << lv;
     double length = 0, width = 0;
     // Set length to be the largest size, width the smallest
     std::sort(paras.begin(), paras.end());
     length = 2.0 * paras.back();
     width = 2.0 * paras.front();
     xtalLMap.insert(std::pair<G4LogicalVolume *, Doubles>(lv, Doubles(length, width)));
     dodet = fv.next();
   }
   LogDebug("EcalSim") << "DreamSD: Length Table for " << attribute << " = " << sd << ":";
   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();
     LogDebug("EcalSim") << " " << i << " " << ite->first << " " << name << " L = " << ite->second.first
                         << " W = " << ite->second.second;
   }
 }

void DreamSD::initRun ( )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 74 of file DreamSD.cc.

References materialPropertiesTable, setPbWO2MaterialProperties_(), and xtalLMap.

Referenced by ~DreamSD().

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

uint32_t DreamSD::setDetUnitId ( const G4Step * aStep )

overridevirtual

Implements CaloSD.

Definition at line 92 of file DreamSD.cc.

References triggerObjects_cff::id, LogDebug, readBothSide_, and side.

Referenced by ~DreamSD().

                                                   {
   const G4VTouchable *touch = aStep->GetPreStepPoint()->GetTouchable();
   uint32_t id = (touch->GetReplicaNumber(1)) * 10 + (touch->GetReplicaNumber(0));
   side = readBothSide_ ? -1 : 1;
   if (side < 0) {
     ++id;
   }
   LogDebug("EcalSim") << "DreamSD:: ID " << id;
   return id;
 }

bool DreamSD::setPbWO2MaterialProperties_ ( G4Material * aMaterial )

private

Sets material properties at run-time...

Definition at line 371 of file DreamSD.cc.

References chAngleIntegrals_, LogDebug, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by initRun().

                                                                {
   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_.reset(new 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;
   }
 
   LogDebug("EcalSim") << "Material properties set for " << aMaterial->GetName();
 
   return true;
 }