CMS 3D CMS Logo

List of all members | Public Member Functions | Protected Member Functions | Private Types | Private Member Functions | Private Attributes | Static Private Attributes
DreamSD Class Reference

#include <DreamSD.h>

Inheritance diagram for DreamSD:
CaloSD SensitiveCaloDetector Observer< const BeginOfRun *> Observer< const BeginOfEvent *> Observer< const BeginOfTrack *> Observer< const EndOfTrack *> Observer< const EndOfEvent *> SensitiveDetector

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)
 
CaloG4HitcreateNewHit (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
TrackInformationcmsTrackInformation (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 *, DoublesDimensionMap
 
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::DDCompactViewcpvDD4hep_
 
const DDCompactViewcpvDDD_
 
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]
 
CaloG4HitcurrentHit [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(), dqmdumpme::k, Skims_PA_cff::name, AlCaHLTBitMon_ParallelJobs::p, readBothSide_, slopeLY_, and useBirk_.

38  : CaloSD(name, clg, p, manager), cpvDDD_(cpvDDD), cpvDD4hep_(cpvDD4hep) {
39  edm::ParameterSet m_EC = p.getParameter<edm::ParameterSet>("ECalSD");
40  useBirk_ = m_EC.getParameter<bool>("UseBirkLaw");
41  doCherenkov_ = m_EC.getParameter<bool>("doCherenkov");
42  birk1_ = m_EC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
43  birk2_ = m_EC.getParameter<double>("BirkC2");
44  birk3_ = m_EC.getParameter<double>("BirkC3");
45  slopeLY_ = m_EC.getParameter<double>("SlopeLightYield");
46  readBothSide_ = m_EC.getUntrackedParameter<bool>("ReadBothSide", false);
47  dd4hep_ = p.getParameter<bool>("g4GeometryDD4hepSource");
48 
49  chAngleIntegrals_.reset(nullptr);
50 
51  edm::LogVerbatim("EcalSim") << "Constructing a DreamSD with name " << GetName()
52  << "\nDreamSD:: Use of Birks law is set to " << useBirk_
53  << " with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_
54  << "\n Slope for Light yield is set to " << slopeLY_
55  << "\n Parameterization of Cherenkov is set to " << doCherenkov_
56  << ", readout both sides is " << readBothSide_ << " and dd4hep flag " << dd4hep_;
57 #ifdef EDM_ML_DEBUG
58  edm::LogVerbatim("EcalSim") << GetName() << " initialized";
59  const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
60  unsigned int k(0);
61  for (auto lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite, ++k)
62  edm::LogVerbatim("EcalSim") << "Volume[" << k << "] " << (*lvcite)->GetName();
63 #endif
64  initMap(name);
65 }
const DDCompactView * cpvDDD_
Definition: DreamSD.h:61
Log< level::Info, true > LogVerbatim
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
bool readBothSide_
Definition: DreamSD.h:64
double birk3_
Definition: DreamSD.h:65
bool useBirk_
Definition: DreamSD.h:64
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
T getUntrackedParameter(std::string const &, T const &) const
bool dd4hep_
Definition: DreamSD.h:64
std::unique_ptr< G4PhysicsFreeVector > chAngleIntegrals_
Table of Cherenkov angle integrals vs photon momentum.
Definition: DreamSD.h:72
double birk2_
Definition: DreamSD.h:65
CaloSD(const std::string &aSDname, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, float timeSlice=1., bool ignoreTkID=false, const std::string &newcolname="")
Definition: CaloSD.cc:33
const cms::DDCompactView * cpvDD4hep_
Definition: DreamSD.h:62
bool doCherenkov_
Definition: DreamSD.h:64
double birk1_
Definition: DreamSD.h:65
double slopeLY_
Definition: DreamSD.h:66
void initMap(const std::string &)
Definition: DreamSD.cc:121

◆ ~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_2024v12_cff::beta, ALCARECOTkAlJpsiMuMu_cff::charge, funct::cos(), l1tSlwPFPuppiJets_cfi::cosPhi, Calorimetry_cff::dp, getAverageNumberOfPhotons_(), getPhotonEnergyDeposit_(), materialPropertiesTable_, phi, ElectronMcFakeValidator_cfi::Pmax, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, funct::sin(), l1tSlwPFPuppiJets_cfi::sinPhi, and mathSSE::sqrt().

Referenced by getEnergyDeposit().

243  {
244  double cherenkovEnergy = 0;
246  return cherenkovEnergy;
247  G4Material *material = aStep->GetTrack()->GetMaterial();
248 
249  // Retrieve refractive index
250  G4MaterialPropertyVector *Rindex = materialPropertiesTable_->GetProperty("RINDEX");
251  if (Rindex == nullptr) {
252  edm::LogWarning("EcalSim") << "Couldn't retrieve refractive index";
253  return cherenkovEnergy;
254  }
255 
256  // V.Ivanchenko - temporary close log output for 9.5
257  // Material refraction properties
258  int Rlength = Rindex->GetVectorLength() - 1;
259  double Pmin = Rindex->Energy(0);
260  double Pmax = Rindex->Energy(Rlength);
261 #ifdef EDM_ML_DEBUG
262  edm::LogVerbatim("EcalSim") << "Material properties: \n Pmin = " << Pmin << " Pmax = " << Pmax;
263 #endif
264  // Get particle properties
265  const G4StepPoint *pPreStepPoint = aStep->GetPreStepPoint();
266  const G4StepPoint *pPostStepPoint = aStep->GetPostStepPoint();
267  const G4ThreeVector &x0 = pPreStepPoint->GetPosition();
268  G4ThreeVector p0 = aStep->GetDeltaPosition().unit();
269  const G4DynamicParticle *aParticle = aStep->GetTrack()->GetDynamicParticle();
270  const double charge = aParticle->GetDefinition()->GetPDGCharge();
271  // beta is averaged over step
272  double beta = 0.5 * (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta());
273  double BetaInverse = 1.0 / beta;
274 
275 #ifdef EDM_ML_DEBUG
276  edm::LogVerbatim("EcalSim") << "Particle properties: \n charge = " << charge << " beta = " << beta;
277 #endif
278 
279  // Now get number of photons generated in this step
280  double meanNumberOfPhotons = getAverageNumberOfPhotons_(charge, beta, material, Rindex);
281  if (meanNumberOfPhotons <= 0.0) { // Don't do anything
282 #ifdef EDM_ML_DEBUG
283  edm::LogVerbatim("EcalSim") << "Mean number of photons is zero: " << meanNumberOfPhotons << ", stopping here";
284 #endif
285  return cherenkovEnergy;
286  }
287 
288  // number of photons is in unit of Geant4...
289  meanNumberOfPhotons *= aStep->GetStepLength();
290 
291  // Now get a poisson distribution
292  int numPhotons = static_cast<int>(G4Poisson(meanNumberOfPhotons));
293  // edm::LogVerbatim("EcalSim") << "Number of photons = " << numPhotons;
294  if (numPhotons <= 0) {
295 #ifdef EDM_ML_DEBUG
296  edm::LogVerbatim("EcalSim") << "Poission number of photons is zero: " << numPhotons << ", stopping here";
297 #endif
298  return cherenkovEnergy;
299  }
300 
301  // Material refraction properties
302  double dp = Pmax - Pmin;
303  double maxCos = BetaInverse / (*Rindex)[Rlength];
304  double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
305 
306  // Finally: get contribution of each photon
307  for (int iPhoton = 0; iPhoton < numPhotons; ++iPhoton) {
308  // Determine photon momentum
309  double randomNumber;
310  double sampledMomentum, sampledRI;
311  double cosTheta, sin2Theta;
312 
313  // sample a momentum (not sure why this is needed!)
314  do {
315  randomNumber = G4UniformRand();
316  sampledMomentum = Pmin + randomNumber * dp;
317  sampledRI = Rindex->Value(sampledMomentum);
318  cosTheta = BetaInverse / sampledRI;
319 
320  sin2Theta = (1.0 - cosTheta) * (1.0 + cosTheta);
321  randomNumber = G4UniformRand();
322 
323  } while (randomNumber * maxSin2 > sin2Theta);
324 
325  // Generate random position of photon on cone surface
326  // defined by Theta
327  randomNumber = G4UniformRand();
328 
329  double phi = twopi * randomNumber;
330  double sinPhi = sin(phi);
331  double cosPhi = cos(phi);
332 
333  // Create photon momentum direction vector
334  // The momentum direction is still w.r.t. the coordinate system where the
335  // primary particle direction is aligned with the z axis
336  double sinTheta = sqrt(sin2Theta);
337  double px = sinTheta * cosPhi;
338  double py = sinTheta * sinPhi;
339  double pz = cosTheta;
340  G4ThreeVector photonDirection(px, py, pz);
341 
342  // Rotate momentum direction back to global (crystal) reference system
343  photonDirection.rotateUz(p0);
344 
345  // Create photon position and momentum
346  randomNumber = G4UniformRand();
347  G4ThreeVector photonPosition = x0 + randomNumber * aStep->GetDeltaPosition();
348  G4ThreeVector photonMomentum = sampledMomentum * photonDirection;
349 
350  // Collect energy on APD
351  cherenkovEnergy += getPhotonEnergyDeposit_(photonMomentum, photonPosition, aStep);
352  }
353  return cherenkovEnergy;
354 }
Log< level::Info, true > LogVerbatim
double getPhotonEnergyDeposit_(const G4ParticleMomentum &p, const G4ThreeVector &x, const G4Step *aStep)
Returns energy deposit for a given photon.
Definition: DreamSD.cc:504
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
T sqrt(T t)
Definition: SSEVec.h:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
double getAverageNumberOfPhotons_(const double charge, const double beta, const G4Material *aMaterial, const G4MaterialPropertyVector *rIndex)
Returns average number of photons created by track.
Definition: DreamSD.cc:359
G4MaterialPropertiesTable * materialPropertiesTable_
Definition: DreamSD.h:73
Log< level::Warning, false > LogWarning

◆ crystalLength()

double DreamSD::crystalLength ( G4LogicalVolume *  lv) const
private

Definition at line 223 of file DreamSD.cc.

References xtalLMap_.

Referenced by curve_LY(), and getPhotonEnergyDeposit_().

223  {
224  double length = -1.;
225  DimensionMap::const_iterator ite = xtalLMap_.find(lv);
226  if (ite != xtalLMap_.end())
227  length = ite->second.first;
228  return length;
229 }
DimensionMap xtalLMap_
Definition: DreamSD.h:67

◆ 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_().

232  {
233  double width = -1.;
234  DimensionMap::const_iterator ite = xtalLMap_.find(lv);
235  if (ite != xtalLMap_.end())
236  width = ite->second.second;
237  return width;
238 }
DimensionMap xtalLMap_
Definition: DreamSD.h:67

◆ 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().

196  {
197  auto const stepPoint = aStep->GetPreStepPoint();
198  auto const lv = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
199  G4String nameVolume = lv->GetName();
200 
201  double weight = 1.;
202  G4ThreeVector localPoint = setToLocal(stepPoint->GetPosition(), stepPoint->GetTouchable());
203  double crlength = crystalLength(lv);
204  double localz = localPoint.x();
205  double dapd = 0.5 * crlength - flag * localz; // Distance from closest APD
206  if (dapd >= -0.1 || dapd <= crlength + 0.1) {
207  if (dapd <= 100.)
208  weight = 1.0 + slopeLY_ - dapd * 0.01 * slopeLY_;
209  } else {
210  edm::LogWarning("EcalSim") << "DreamSD: light coll curve : wrong distance "
211  << "to APD " << dapd << " crlength = " << crlength << " crystal name = " << nameVolume
212  << " z of localPoint = " << localz << " take weight = " << weight;
213  }
214 #ifdef EDM_ML_DEBUG
215  edm::LogVerbatim("EcalSim") << "DreamSD, light coll curve : " << dapd << " crlength = " << crlength
216  << " crystal name = " << nameVolume << " z of localPoint = " << localz
217  << " take weight = " << weight;
218 #endif
219  return weight;
220 }
Log< level::Info, true > LogVerbatim
double crystalLength(G4LogicalVolume *) const
Definition: DreamSD.cc:223
G4ThreeVector setToLocal(const G4ThreeVector &, const G4VTouchable *) const
Definition: CaloSD.cc:454
Definition: weight.py:1
double slopeLY_
Definition: DreamSD.h:66
Log< level::Warning, false > LogWarning

◆ fillMap()

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

Definition at line 176 of file DreamSD.cc.

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

Referenced by initMap().

176  {
177  const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
178  edm::LogVerbatim("EcalSim") << "LV Store with " << lvs->size() << " elements";
179  G4LogicalVolume *lv = nullptr;
180  for (auto lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
181  edm::LogVerbatim("EcalSim") << name << " vs " << (*lvcite)->GetName();
182  std::string namex = static_cast<std::string>((*lvcite)->GetName());
183  if (name == DD4hep2DDDName::noNameSpace(static_cast<std::string>(namex))) {
184  lv = (*lvcite);
185  break;
186  }
187  }
188 #ifdef EDM_ML_DEBUG
189  edm::LogVerbatim("EcalSim") << "DreamSD::fillMap (for " << name << " Logical Volume " << lv << " Length " << length
190  << " Width " << width;
191 #endif
192  xtalLMap_.insert(std::pair<G4LogicalVolume *, Doubles>(lv, Doubles(length, width)));
193 }
Log< level::Info, true > LogVerbatim
DimensionMap xtalLMap_
Definition: DreamSD.h:67
std::pair< double, double > Doubles
Definition: DreamSD.h:37
std::string noNameSpace(const std::string &name)

◆ 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_2024v12_cff::beta, chAngleIntegrals_, ALCARECOTkAlJpsiMuMu_cff::charge, Calorimetry_cff::dp, metDiagnosticParameterSet_cfi::nMax, metDiagnosticParameterSet_cfi::nMin, and ElectronMcFakeValidator_cfi::Pmax.

Referenced by cherenkovDeposit_().

362  {
363  const double rFact = 369.81 / (eV * cm);
364 
365  if (beta <= 0.0)
366  return 0.0;
367 
368  double BetaInverse = 1. / beta;
369 
370  // Vectors used in computation of Cerenkov Angle Integral:
371  // - Refraction Indices for the current material
372  // - new G4PhysicsOrderedFreeVector allocated to hold CAI's
373 
374  // Min and Max photon momenta
375  int Rlength = Rindex->GetVectorLength() - 1;
376  double Pmin = Rindex->Energy(0);
377  double Pmax = Rindex->Energy(Rlength);
378 
379  // Min and Max Refraction Indices
380  double nMin = (*Rindex)[0];
381  double nMax = (*Rindex)[Rlength];
382 
383  // Max Cerenkov Angle Integral
384  double CAImax = chAngleIntegrals_.get()->GetMaxEnergy();
385 
386  double dp = 0., ge = 0., CAImin = 0.;
387 
388  // If n(Pmax) < 1/Beta -- no photons generated
389  if (nMax < BetaInverse) {
390  }
391 
392  // otherwise if n(Pmin) >= 1/Beta -- photons generated
393  else if (nMin > BetaInverse) {
394  dp = Pmax - Pmin;
395  ge = CAImax;
396  }
397  // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
398  // we need to find a P such that the value of n(P) == 1/Beta.
399  // Interpolation is performed by the GetPhotonEnergy() and
400  // GetProperty() methods of the G4MaterialPropertiesTable and
401  // the GetValue() method of G4PhysicsVector.
402  else {
403  Pmin = Rindex->Value(BetaInverse);
404  dp = Pmax - Pmin;
405  // need boolean for current implementation of G4PhysicsVector
406  // ==> being phased out
407  double CAImin = chAngleIntegrals_->Value(Pmin);
408  ge = CAImax - CAImin;
409  }
410 
411  // Calculate number of photons
412  double numPhotons = rFact * charge / eplus * charge / eplus * (dp - ge * BetaInverse * BetaInverse);
413 
414  edm::LogVerbatim("EcalSim") << "@SUB=getAverageNumberOfPhotons\nCAImin = " << CAImin << "\nCAImax = " << CAImax
415  << "\ndp = " << dp << ", ge = " << ge << "\nnumPhotons = " << numPhotons;
416  return numPhotons;
417 }
Log< level::Info, true > LogVerbatim
std::unique_ptr< G4PhysicsFreeVector > chAngleIntegrals_
Table of Cherenkov angle integrals vs photon momentum.
Definition: DreamSD.h:72

◆ 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.

68  {
69  // take into account light collection curve for crystals
70  double weight = curve_LY(aStep, side_);
71  if (useBirk_)
73  double edep = aStep->GetTotalEnergyDeposit() * weight;
74 
75  // Get Cerenkov contribution
76  if (doCherenkov_) {
77  edep += cherenkovDeposit_(aStep);
78  }
79 #ifdef EDM_ML_DEBUG
80  edm::LogVerbatim("EcalSim") << "DreamSD:: " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() << " Side "
81  << side_ << " Light Collection Efficiency " << weight << " Weighted Energy Deposit "
82  << edep / CLHEP::MeV << " MeV";
83 #endif
84  return edep;
85 }
Log< level::Info, true > LogVerbatim
Definition: weight.py:1
double birk3_
Definition: DreamSD.h:65
bool useBirk_
Definition: DreamSD.h:64
double cherenkovDeposit_(const G4Step *aStep)
Returns the total energy due to Cherenkov radiation.
Definition: DreamSD.cc:243
int side_
Definition: DreamSD.h:69
double curve_LY(const G4Step *, int)
Definition: DreamSD.cc:196
double birk2_
Definition: DreamSD.h:65
double getAttenuation(const G4Step *aStep, double birk1, double birk2, double birk3) const
Definition: CaloSD.cc:727
bool doCherenkov_
Definition: DreamSD.h:64
double birk1_
Definition: DreamSD.h:65

◆ getPhotonEnergyDeposit_()

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

Returns energy deposit for a given photon.

Definition at line 504 of file DreamSD.cc.

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

Referenced by cherenkovDeposit_().

504  {
505  double energy = 0;
506 
507  // Crystal dimensions
508 
509  // edm::LogVerbatim("EcalSim") << p << x;
510 
511  // 1. Check if this photon goes straight to the APD:
512  // - assume that APD is at x=xtalLength/2.0
513  // - extrapolate from x=x0 to x=xtalLength/2.0 using momentum in x-y
514 
515  G4StepPoint *stepPoint = aStep->GetPreStepPoint();
516  G4LogicalVolume *lv = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
517  G4String nameVolume = lv->GetName();
518 
519  double crlength = crystalLength(lv);
520  double crwidth = crystalWidth(lv);
521  double dapd = 0.5 * crlength - x.x(); // Distance from closest APD
522  double y = p.y() / p.x() * dapd;
523 
524 #ifdef EDM_ML_DEBUG
525  edm::LogVerbatim("EcalSim") << "Distance to APD: " << dapd << " - y at APD: " << y;
526 #endif
527  // Not straight: compute probability
528  if (std::abs(y) > crwidth * 0.5) {
529  }
530 
531  // 2. Retrieve efficiency for this wavelength (in nm, from MeV)
532  double waveLength = p.mag() * 1.239e8;
533 
534  energy = p.mag() * PMTResponse::getEfficiency(waveLength);
535 #ifdef EDM_ML_DEBUG
536  edm::LogVerbatim("EcalSim") << "Wavelength: " << waveLength << " - Energy: " << energy;
537 #endif
538  return energy;
539 }
Log< level::Info, true > LogVerbatim
double crystalLength(G4LogicalVolume *) const
Definition: DreamSD.cc:223
double crystalWidth(G4LogicalVolume *) const
Definition: DreamSD.cc:232
static double getEfficiency(const double &waveLengthNm)
Return efficiency for given photon wavelength (in nm)
Definition: PMTResponse.cc:6
Abs< T >::type abs(const T &t)
Definition: Abs.h:22

◆ 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(), Skims_PA_cff::name, cms::dd::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().

121  {
122  if (dd4hep_) {
123  const cms::DDFilter filter("ReadOutName", sd);
125  while (fv.firstChild()) {
126  std::string name = DD4hep2DDDName::noNameSpace(static_cast<std::string>(fv.name()));
127  std::vector<double> paras(fv.parameters());
128 #ifdef EDM_ML_DEBUG
129  edm::LogVerbatim("EcalSim") << "DreamSD::initMap (for " << sd << "): Solid " << name << " Shape "
130  << cms::dd::name(cms::DDSolidShapeMap, fv.shape()) << " Parameter 0 = " << paras[0];
131 #endif
132  // Set length to be the largest size, width the smallest
133  std::sort(paras.begin(), paras.end());
134  double length = 2.0 * k_ScaleFromDD4hepToG4 * paras.back();
135  double width = 2.0 * k_ScaleFromDD4hepToG4 * paras.front();
136  fillMap(name, length, width);
137  }
138  } else {
139  DDSpecificsMatchesValueFilter filter{DDValue("ReadOutName", sd, 0)};
140  DDFilteredView fv((*cpvDDD_), filter);
141  fv.firstChild();
142  bool dodet = true;
143  while (dodet) {
144  const DDSolid &sol = fv.logicalPart().solid();
145  std::vector<double> paras(sol.parameters());
146  std::string name = static_cast<std::string>(sol.name().name());
147 #ifdef EDM_ML_DEBUG
148  edm::LogVerbatim("EcalSim") << "DreamSD::initMap (for " << sd << "): Solid " << name << " Shape " << sol.shape()
149  << " Parameter 0 = " << paras[0];
150 #endif
151  // Set length to be the largest size, width the smallest
152  std::sort(paras.begin(), paras.end());
153  double length = 2.0 * k_ScaleFromDDDToG4 * paras.back();
154  double width = 2.0 * k_ScaleFromDDDToG4 * paras.front();
155  fillMap(name, length, width);
156  dodet = fv.next();
157  }
158  }
159 #ifdef EDM_ML_DEBUG
160  edm::LogVerbatim("EcalSim") << "DreamSD: Length Table for ReadOutName = " << sd << ":";
161 #endif
162  DimensionMap::const_iterator ite = xtalLMap_.begin();
163  int i = 0;
164  for (; ite != xtalLMap_.end(); ite++, i++) {
165  G4String name = "Unknown";
166  if (ite->first != nullptr)
167  name = (ite->first)->GetName();
168 #ifdef EDM_ML_DEBUG
169  edm::LogVerbatim("EcalSim") << " " << i << " " << ite->first << " " << name << " L = " << ite->second.first
170  << " W = " << ite->second.second;
171 #endif
172  }
173 }
const DDCompactView * cpvDDD_
Definition: DreamSD.h:61
Log< level::Info, true > LogVerbatim
static constexpr double k_ScaleFromDD4hepToG4
Definition: DreamSD.h:59
DimensionMap xtalLMap_
Definition: DreamSD.h:67
std::string name(Mapping a, V value)
Definition: DDSolidShapes.h:31
void fillMap(const std::string &, double, double)
Definition: DreamSD.cc:176
static constexpr double k_ScaleFromDDDToG4
Definition: DreamSD.h:58
A DDSolid represents the shape of a part.
Definition: DDSolid.h:39
bool dd4hep_
Definition: DreamSD.h:64
const std::array< const cms::dd::NameValuePair< DDSolidShape >, 21 > DDSolidShapeMap
Definition: DDSolidShapes.h:99
constexpr float sol
Definition: Config.h:13
const cms::DDCompactView * cpvDD4hep_
Definition: DreamSD.h:62
std::string noNameSpace(const std::string &name)

◆ initRun()

void DreamSD::initRun ( )
overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 88 of file DreamSD.cc.

References materialPropertiesTable_, setPbWO2MaterialProperties_(), and xtalLMap_.

88  {
89  // Get the material and set properties if needed
90  DimensionMap::const_iterator ite = xtalLMap_.begin();
91  const G4LogicalVolume *lv = (ite->first);
92  G4Material *material = lv->GetMaterial();
93 #ifdef EDM_ML_DEBUG
94  edm::LogVerbatim("EcalSim") << "DreamSD::initRun: Initializes for material " << material->GetName() << " in "
95  << lv->GetName();
96 #endif
97  materialPropertiesTable_ = material->GetMaterialPropertiesTable();
99  if (!setPbWO2MaterialProperties_(material)) {
100  edm::LogWarning("EcalSim") << "Couldn't retrieve material properties table\n Material = " << material->GetName();
101  }
102  materialPropertiesTable_ = material->GetMaterialPropertiesTable();
103  }
104 }
Log< level::Info, true > LogVerbatim
DimensionMap xtalLMap_
Definition: DreamSD.h:67
bool setPbWO2MaterialProperties_(G4Material *aMaterial)
Sets material properties at run-time...
Definition: DreamSD.cc:422
G4MaterialPropertiesTable * materialPropertiesTable_
Definition: DreamSD.h:73
Log< level::Warning, false > LogWarning

◆ setDetUnitId()

uint32_t DreamSD::setDetUnitId ( const G4Step *  aStep)
overridevirtual

Implements CaloSD.

Definition at line 107 of file DreamSD.cc.

References l1ctLayer2EG_cff::id, readBothSide_, and side_.

107  {
108  const G4VTouchable *touch = aStep->GetPreStepPoint()->GetTouchable();
109  uint32_t id = (touch->GetReplicaNumber(1)) * 10 + (touch->GetReplicaNumber(0));
110  side_ = readBothSide_ ? -1 : 1;
111  if (side_ < 0) {
112  ++id;
113  }
114 #ifdef EDM_ML_DEBUG
115  edm::LogVerbatim("EcalSim") << "DreamSD:: ID " << id;
116 #endif
117  return id;
118 }
Log< level::Info, true > LogVerbatim
bool readBothSide_
Definition: DreamSD.h:64
int side_
Definition: DreamSD.h:69

◆ setPbWO2MaterialProperties_()

bool DreamSD::setPbWO2MaterialProperties_ ( G4Material *  aMaterial)
private

Sets material properties at run-time...

Definition at line 422 of file DreamSD.cc.

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

Referenced by initRun().

422  {
423  std::string pbWO2Name("E_PbWO4");
424  std::string name = static_cast<std::string>(aMaterial->GetName());
425  if (DD4hep2DDDName::noNameSpace(name) != pbWO2Name) { // Wrong material!
426  edm::LogWarning("EcalSim") << "This is not the right material: "
427  << "expecting " << pbWO2Name << ", got " << aMaterial->GetName();
428  return false;
429  }
430 
431  G4MaterialPropertiesTable *table = new G4MaterialPropertiesTable();
432 
433  // Refractive index as a function of photon momentum
434  // FIXME: Should somehow put that in the configuration
435  const int nEntries = 14;
436  double PhotonEnergy[nEntries] = {1.7712 * eV,
437  1.8368 * eV,
438  1.90745 * eV,
439  1.98375 * eV,
440  2.0664 * eV,
441  2.15625 * eV,
442  2.25426 * eV,
443  2.3616 * eV,
444  2.47968 * eV,
445  2.61019 * eV,
446  2.75521 * eV,
447  2.91728 * eV,
448  3.09961 * eV,
449  3.30625 * eV};
450  double RefractiveIndex[nEntries] = {2.17728,
451  2.18025,
452  2.18357,
453  2.18753,
454  2.19285,
455  2.19813,
456  2.20441,
457  2.21337,
458  2.22328,
459  2.23619,
460  2.25203,
461  2.27381,
462  2.30282,
463  2.34666};
464 
465  table->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex, nEntries);
466  aMaterial->SetMaterialPropertiesTable(table); // FIXME: could this leak? What does G4 do?
467 
468  // Calculate Cherenkov angle integrals:
469  // This is an ad-hoc solution (we hold it in the class, not in the material)
470  chAngleIntegrals_ = std::make_unique<G4PhysicsFreeVector>(nEntries);
471 
472  int index = 0;
473  double currentRI = RefractiveIndex[index];
474  double currentPM = PhotonEnergy[index];
475  double currentCAI = 0.0;
476  chAngleIntegrals_.get()->PutValue(0, currentPM, currentCAI);
477  double prevPM = currentPM;
478  double prevCAI = currentCAI;
479  double prevRI = currentRI;
480  while (++index < nEntries) {
481  currentRI = RefractiveIndex[index];
482  currentPM = PhotonEnergy[index];
483  currentCAI = 0.5 * (1.0 / (prevRI * prevRI) + 1.0 / (currentRI * currentRI));
484  currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
485 
486  chAngleIntegrals_.get()->PutValue(index, currentPM, currentCAI);
487 
488  prevPM = currentPM;
489  prevCAI = currentCAI;
490  prevRI = currentRI;
491  }
492 
493 #ifdef EDM_ML_DEBUG
494  edm::LogVerbatim("EcalSim") << "Material properties set for " << aMaterial->GetName();
495 #endif
496  return true;
497 }
Log< level::Info, true > LogVerbatim
std::unique_ptr< G4PhysicsFreeVector > chAngleIntegrals_
Table of Cherenkov angle integrals vs photon momentum.
Definition: DreamSD.h:72
std::string noNameSpace(const std::string &name)
Log< level::Warning, false > LogWarning

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().