CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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)
 
void clear () override
 
void clearHits () override
 
void DrawAll () override
 
void EndOfEvent (G4HCofThisEvent *eventHC) override
 
void fillHits (edm::PCaloHitContainer &, const std::string &) override
 
void Initialize (G4HCofThisEvent *HCE) override
 
bool isItFineCalo (const G4VTouchable *touch)
 
void PrintAll () override
 
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)
 
- 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)
 
 ~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 ()
 
CaloG4HitcreateNewHit (const G4Step *, const G4Track *)
 
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 ()
 
double getResponseWt (const G4Track *)
 
virtual int getTrackID (const G4Track *)
 
bool hitExists (const G4Step *)
 
void ignoreRejection ()
 
virtual void initEvent (const BeginOfEvent *)
 
void printDetectorLevels (const G4VTouchable *) const
 
void processHit (const G4Step *step)
 
void resetForNewPrimary (const G4Step *)
 
void setNumberCheckedHits (int val)
 
void setParameterized (bool val)
 
G4ThreeVector setToGlobal (const G4ThreeVector &, const G4VTouchable *) const
 
G4ThreeVector setToLocal (const G4ThreeVector &, const G4VTouchable *) const
 
virtual int setTrackID (const G4Step *)
 
void setUseMap (bool val)
 
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 *)
 
- 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 *, 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::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
CaloG4HitcurrentHit
 
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 22 of file DreamSD.h.

Member Typedef Documentation

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

Definition at line 40 of file DreamSD.h.

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

Definition at line 39 of file DreamSD.h.

Constructor & Destructor Documentation

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

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

37  : CaloSD(name, clg, p, manager), cpvDDD_(cpvDDD), cpvDD4hep_(cpvDD4hep) {
38  edm::ParameterSet m_EC = p.getParameter<edm::ParameterSet>("ECalSD");
39  useBirk_ = m_EC.getParameter<bool>("UseBirkLaw");
40  doCherenkov_ = m_EC.getParameter<bool>("doCherenkov");
41  birk1_ = m_EC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
42  birk2_ = m_EC.getParameter<double>("BirkC2");
43  birk3_ = m_EC.getParameter<double>("BirkC3");
44  slopeLY_ = m_EC.getParameter<double>("SlopeLightYield");
45  readBothSide_ = m_EC.getUntrackedParameter<bool>("ReadBothSide", false);
46  dd4hep_ = p.getParameter<bool>("g4GeometryDD4hepSource");
47 
48  chAngleIntegrals_.reset(nullptr);
49 
50  edm::LogVerbatim("EcalSim") << "Constructing a DreamSD with name " << GetName()
51  << "\nDreamSD:: Use of Birks law is set to " << useBirk_
52  << " with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_
53  << "\n Slope for Light yield is set to " << slopeLY_
54  << "\n Parameterization of Cherenkov is set to " << doCherenkov_
55  << ", readout both sides is " << readBothSide_ << " and dd4hep flag " << dd4hep_;
56 #ifdef EDM_ML_DEBUG
57  edm::LogVerbatim("EcalSim") << GetName() << " initialized";
58  const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
59  unsigned int k(0);
60  for (auto lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite, ++k)
61  edm::LogVerbatim("EcalSim") << "Volume[" << k << "] " << (*lvcite)->GetName();
62 #endif
63  initMap(name);
64 }
const DDCompactView * cpvDDD_
Definition: DreamSD.h:63
Log< level::Info, true > LogVerbatim
T getUntrackedParameter(std::string const &, T const &) const
bool readBothSide_
Definition: DreamSD.h:66
double birk3_
Definition: DreamSD.h:67
bool useBirk_
Definition: DreamSD.h:66
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
bool dd4hep_
Definition: DreamSD.h:66
std::unique_ptr< G4PhysicsFreeVector > chAngleIntegrals_
Table of Cherenkov angle integrals vs photon momentum.
Definition: DreamSD.h:74
double birk2_
Definition: DreamSD.h:67
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
const cms::DDCompactView * cpvDD4hep_
Definition: DreamSD.h:64
bool doCherenkov_
Definition: DreamSD.h:66
CaloSD(const std::string &aSDname, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, float timeSlice=1., bool ignoreTkID=false)
Definition: CaloSD.cc:33
double birk1_
Definition: DreamSD.h:67
double slopeLY_
Definition: DreamSD.h:68
void initMap(const std::string &)
Definition: DreamSD.cc:120
DreamSD::~DreamSD ( )
inlineoverride

Definition at line 30 of file DreamSD.h.

30 {}

Member Function Documentation

double DreamSD::cherenkovDeposit_ ( const G4Step *  aStep)
private

Returns the total energy due to Cherenkov radiation.

Definition at line 242 of file DreamSD.cc.

References HLT_FULL_cff::beta, RecoTauCleanerPlugins::charge, funct::cos(), Phase1L1TJetProducer_cfi::cosPhi, getAverageNumberOfPhotons_(), getPhotonEnergyDeposit_(), materialPropertiesTable_, phi, funct::sin(), Phase1L1TJetProducer_cfi::sinPhi, and mathSSE::sqrt().

Referenced by getEnergyDeposit().

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

Definition at line 222 of file DreamSD.cc.

References xtalLMap_.

Referenced by curve_LY(), and getPhotonEnergyDeposit_().

222  {
223  double length = -1.;
224  DimensionMap::const_iterator ite = xtalLMap_.find(lv);
225  if (ite != xtalLMap_.end())
226  length = ite->second.first;
227  return length;
228 }
DimensionMap xtalLMap_
Definition: DreamSD.h:69
double DreamSD::crystalWidth ( G4LogicalVolume *  lv) const
private

Definition at line 231 of file DreamSD.cc.

References xtalLMap_.

Referenced by getPhotonEnergyDeposit_().

231  {
232  double width = -1.;
233  DimensionMap::const_iterator ite = xtalLMap_.find(lv);
234  if (ite != xtalLMap_.end())
235  width = ite->second.second;
236  return width;
237 }
DimensionMap xtalLMap_
Definition: DreamSD.h:69
double DreamSD::curve_LY ( const G4Step *  aStep,
int  flag 
)
private

Definition at line 195 of file DreamSD.cc.

References crystalLength(), CaloSD::setToLocal(), slopeLY_, and histoStyle::weight.

Referenced by getEnergyDeposit().

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

Definition at line 175 of file DreamSD.cc.

References AlCaHLTBitMon_QueryRunRegistry::string, and xtalLMap_.

Referenced by initMap().

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

References HLT_FULL_cff::beta, and chAngleIntegrals_.

Referenced by cherenkovDeposit_().

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

Reimplemented from CaloSD.

Definition at line 67 of file DreamSD.cc.

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

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

Returns energy deposit for a given photon.

Definition at line 503 of file DreamSD.cc.

References funct::abs(), crystalLength(), crystalWidth(), relval_parameters_module::energy, PMTResponse::getEfficiency(), and y.

Referenced by cherenkovDeposit_().

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

Definition at line 120 of file DreamSD.cc.

References cpvDD4hep_, cpvDDD_, dd4hep_, cms::DDSolidShapeMap, fillMap(), alcazmumu_cfi::filter, DDFilteredView::firstChild(), cms::DDFilteredView::firstChild(), mps_fire::i, k_ScaleFromDD4hepToG4, k_ScaleFromDDDToG4, DDFilteredView::logicalPart(), cms::dd::name(), DDName::name(), DDBase< N, C >::name(), mergeVDriftHistosByStation::name, cms::DDFilteredView::name(), DDFilteredView::next(), DDSolid::parameters(), cms::DDFilteredView::parameters(), DDSolid::shape(), cms::DDFilteredView::shape(), mkfit::Const::sol, DDLogicalPart::solid(), AlCaHLTBitMon_QueryRunRegistry::string, and xtalLMap_.

Referenced by DreamSD().

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

Reimplemented from CaloSD.

Definition at line 87 of file DreamSD.cc.

References materialPropertiesTable_, setPbWO2MaterialProperties_(), and xtalLMap_.

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

Implements CaloSD.

Definition at line 106 of file DreamSD.cc.

References gpuClustering::id, readBothSide_, and side_.

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

Sets material properties at run-time...

Definition at line 421 of file DreamSD.cc.

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

Referenced by initRun().

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

Member Data Documentation

double DreamSD::birk1_
private

Definition at line 67 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

double DreamSD::birk2_
private

Definition at line 67 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

double DreamSD::birk3_
private

Definition at line 67 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

std::unique_ptr<G4PhysicsFreeVector> DreamSD::chAngleIntegrals_
private

Table of Cherenkov angle integrals vs photon momentum.

Definition at line 74 of file DreamSD.h.

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

const cms::DDCompactView* DreamSD::cpvDD4hep_
private

Definition at line 64 of file DreamSD.h.

Referenced by initMap().

const DDCompactView* DreamSD::cpvDDD_
private

Definition at line 63 of file DreamSD.h.

Referenced by initMap().

bool DreamSD::dd4hep_
private

Definition at line 66 of file DreamSD.h.

Referenced by DreamSD(), and initMap().

bool DreamSD::doCherenkov_
private

Definition at line 66 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

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

Definition at line 61 of file DreamSD.h.

Referenced by initMap().

constexpr double DreamSD::k_ScaleFromDDDToG4 = 1.0
staticprivate

Definition at line 60 of file DreamSD.h.

Referenced by initMap().

G4MaterialPropertiesTable* DreamSD::materialPropertiesTable_
private

Definition at line 75 of file DreamSD.h.

Referenced by cherenkovDeposit_(), and initRun().

int DreamSD::nphotons_
private

Definition at line 77 of file DreamSD.h.

bool DreamSD::readBothSide_
private

Definition at line 66 of file DreamSD.h.

Referenced by DreamSD(), and setDetUnitId().

int DreamSD::side_
private

Definition at line 71 of file DreamSD.h.

Referenced by getEnergyDeposit(), and setDetUnitId().

double DreamSD::slopeLY_
private

Definition at line 68 of file DreamSD.h.

Referenced by curve_LY(), and DreamSD().

bool DreamSD::useBirk_
private

Definition at line 66 of file DreamSD.h.

Referenced by DreamSD(), and getEnergyDeposit().

DimensionMap DreamSD::xtalLMap_
private

Definition at line 69 of file DreamSD.h.

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