Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-11-09 23:35:30 

0001 #include "DetectorDescription/Core/interface/DDCompactView.h"
0002 #include "DetectorDescription/Core/interface/DDFilter.h"
0003 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0004 #include "DetectorDescription/Core/interface/DDSolid.h"
0005 #include "DetectorDescription/Core/interface/DDSplit.h"
0006 #include "DetectorDescription/Core/interface/DDValue.h"
0007 #include "DetectorDescription/DDCMS/interface/DDCompactView.h"
0008 #include "DetectorDescription/DDCMS/interface/DDFilteredView.h"
0009 #include "SimG4Core/Notification/interface/TrackInformation.h"
0010 
0011 #include "G4LogicalVolume.hh"
0012 #include "G4LogicalVolumeStore.hh"
0013 #include "G4Poisson.hh"
0014 #include "G4Step.hh"
0015 #include "G4Track.hh"
0016 #include "G4VProcess.hh"
0017 
0018 // Histogramming
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 
0021 // Cherenkov
0022 #include "SimG4CMS/CherenkovAnalysis/interface/DreamSD.h"
0023 #include "SimG4CMS/CherenkovAnalysis/interface/PMTResponse.h"
0024 
0025 #include "G4PhysicalConstants.hh"
0026 #include "G4SystemOfUnits.hh"
0027 
0028 //#define EDM_ML_DEBUG
0029 
0030 //________________________________________________________________________________________
0031 DreamSD::DreamSD(const std::string &name,
0032                  const DDCompactView *cpvDDD,
0033                  const cms::DDCompactView *cpvDD4hep,
0034                  const SensitiveDetectorCatalog &clg,
0035                  edm::ParameterSet const &p,
0036                  const SimTrackManager *manager)
0037     : CaloSD(name, clg, p, manager), cpvDDD_(cpvDDD), cpvDD4hep_(cpvDD4hep) {
0038   edm::ParameterSet m_EC = p.getParameter<edm::ParameterSet>("ECalSD");
0039   useBirk_ = m_EC.getParameter<bool>("UseBirkLaw");
0040   doCherenkov_ = m_EC.getParameter<bool>("doCherenkov");
0041   birk1_ = m_EC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
0042   birk2_ = m_EC.getParameter<double>("BirkC2");
0043   birk3_ = m_EC.getParameter<double>("BirkC3");
0044   slopeLY_ = m_EC.getParameter<double>("SlopeLightYield");
0045   readBothSide_ = m_EC.getUntrackedParameter<bool>("ReadBothSide", false);
0046   dd4hep_ = p.getParameter<bool>("g4GeometryDD4hepSource");
0047 
0048   chAngleIntegrals_.reset(nullptr);
0049 
0050   edm::LogVerbatim("EcalSim") << "Constructing a DreamSD  with name " << GetName()
0051                               << "\nDreamSD:: Use of Birks law is set to      " << useBirk_
0052                               << "  with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_
0053                               << "\n          Slope for Light yield is set to " << slopeLY_
0054                               << "\n          Parameterization of Cherenkov is set to " << doCherenkov_
0055                               << ", readout both sides is " << readBothSide_ << " and dd4hep flag " << dd4hep_;
0056 #ifdef EDM_ML_DEBUG
0057   edm::LogVerbatim("EcalSim") << GetName() << " initialized";
0058   const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
0059   unsigned int k(0);
0060   for (auto lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite, ++k)
0061     edm::LogVerbatim("EcalSim") << "Volume[" << k << "] " << (*lvcite)->GetName();
0062 #endif
0063   initMap(name);
0064 }
0065 
0066 //________________________________________________________________________________________
0067 double DreamSD::getEnergyDeposit(const G4Step *aStep) {
0068   // take into account light collection curve for crystals
0069   double weight = curve_LY(aStep, side_);
0070   if (useBirk_)
0071     weight *= getAttenuation(aStep, birk1_, birk2_, birk3_);
0072   double edep = aStep->GetTotalEnergyDeposit() * weight;
0073 
0074   // Get Cerenkov contribution
0075   if (doCherenkov_) {
0076     edep += cherenkovDeposit_(aStep);
0077   }
0078 #ifdef EDM_ML_DEBUG
0079   edm::LogVerbatim("EcalSim") << "DreamSD:: " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() << " Side "
0080                               << side_ << " Light Collection Efficiency " << weight << " Weighted Energy Deposit "
0081                               << edep / CLHEP::MeV << " MeV";
0082 #endif
0083   return edep;
0084 }
0085 
0086 //________________________________________________________________________________________
0087 void DreamSD::initRun() {
0088   // Get the material and set properties if needed
0089   DimensionMap::const_iterator ite = xtalLMap_.begin();
0090   const G4LogicalVolume *lv = (ite->first);
0091   G4Material *material = lv->GetMaterial();
0092 #ifdef EDM_ML_DEBUG
0093   edm::LogVerbatim("EcalSim") << "DreamSD::initRun: Initializes for material " << material->GetName() << " in "
0094                               << lv->GetName();
0095 #endif
0096   materialPropertiesTable_ = material->GetMaterialPropertiesTable();
0097   if (!materialPropertiesTable_) {
0098     if (!setPbWO2MaterialProperties_(material)) {
0099       edm::LogWarning("EcalSim") << "Couldn't retrieve material properties table\n Material = " << material->GetName();
0100     }
0101     materialPropertiesTable_ = material->GetMaterialPropertiesTable();
0102   }
0103 }
0104 
0105 //________________________________________________________________________________________
0106 uint32_t DreamSD::setDetUnitId(const G4Step *aStep) {
0107   const G4VTouchable *touch = aStep->GetPreStepPoint()->GetTouchable();
0108   uint32_t id = (touch->GetReplicaNumber(1)) * 10 + (touch->GetReplicaNumber(0));
0109   side_ = readBothSide_ ? -1 : 1;
0110   if (side_ < 0) {
0111     ++id;
0112   }
0113 #ifdef EDM_ML_DEBUG
0114   edm::LogVerbatim("EcalSim") << "DreamSD:: ID " << id;
0115 #endif
0116   return id;
0117 }
0118 
0119 //________________________________________________________________________________________
0120 void DreamSD::initMap(const std::string &sd) {
0121   if (dd4hep_) {
0122     const cms::DDFilter filter("ReadOutName", sd);
0123     cms::DDFilteredView fv((*cpvDD4hep_), filter);
0124     while (fv.firstChild()) {
0125       std::string name = static_cast<std::string>(dd4hep::dd::noNamespace(fv.name()));
0126       std::vector<double> paras(fv.parameters());
0127 #ifdef EDM_ML_DEBUG
0128       edm::LogVerbatim("EcalSim") << "DreamSD::initMap (for " << sd << "): Solid " << name << " Shape "
0129                                   << cms::dd::name(cms::DDSolidShapeMap, fv.shape()) << " Parameter 0 = " << paras[0];
0130 #endif
0131       // Set length to be the largest size, width the smallest
0132       std::sort(paras.begin(), paras.end());
0133       double length = 2.0 * k_ScaleFromDD4hepToG4 * paras.back();
0134       double width = 2.0 * k_ScaleFromDD4hepToG4 * paras.front();
0135       fillMap(name, length, width);
0136     }
0137   } else {
0138     DDSpecificsMatchesValueFilter filter{DDValue("ReadOutName", sd, 0)};
0139     DDFilteredView fv((*cpvDDD_), filter);
0140     fv.firstChild();
0141     bool dodet = true;
0142     while (dodet) {
0143       const DDSolid &sol = fv.logicalPart().solid();
0144       std::vector<double> paras(sol.parameters());
0145       std::string name = static_cast<std::string>(sol.name().name());
0146 #ifdef EDM_ML_DEBUG
0147       edm::LogVerbatim("EcalSim") << "DreamSD::initMap (for " << sd << "): Solid " << name << " Shape " << sol.shape()
0148                                   << " Parameter 0 = " << paras[0];
0149 #endif
0150       // Set length to be the largest size, width the smallest
0151       std::sort(paras.begin(), paras.end());
0152       double length = 2.0 * k_ScaleFromDDDToG4 * paras.back();
0153       double width = 2.0 * k_ScaleFromDDDToG4 * paras.front();
0154       fillMap(name, length, width);
0155       dodet = fv.next();
0156     }
0157   }
0158 #ifdef EDM_ML_DEBUG
0159   edm::LogVerbatim("EcalSim") << "DreamSD: Length Table for ReadOutName = " << sd << ":";
0160 #endif
0161   DimensionMap::const_iterator ite = xtalLMap_.begin();
0162   int i = 0;
0163   for (; ite != xtalLMap_.end(); ite++, i++) {
0164     G4String name = "Unknown";
0165     if (ite->first != nullptr)
0166       name = (ite->first)->GetName();
0167 #ifdef EDM_ML_DEBUG
0168     edm::LogVerbatim("EcalSim") << " " << i << " " << ite->first << " " << name << " L = " << ite->second.first
0169                                 << " W = " << ite->second.second;
0170 #endif
0171   }
0172 }
0173 
0174 //________________________________________________________________________________________
0175 void DreamSD::fillMap(const std::string &name, double length, double width) {
0176   const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
0177   edm::LogVerbatim("EcalSim") << "LV Store with " << lvs->size() << " elements";
0178   G4LogicalVolume *lv = nullptr;
0179   for (auto lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0180     edm::LogVerbatim("EcalSim") << name << " vs " << (*lvcite)->GetName();
0181     std::string namex = static_cast<std::string>((*lvcite)->GetName());
0182     if (name == static_cast<std::string>(dd4hep::dd::noNamespace(namex))) {
0183       lv = (*lvcite);
0184       break;
0185     }
0186   }
0187 #ifdef EDM_ML_DEBUG
0188   edm::LogVerbatim("EcalSim") << "DreamSD::fillMap (for " << name << " Logical Volume " << lv << " Length " << length
0189                               << " Width " << width;
0190 #endif
0191   xtalLMap_.insert(std::pair<G4LogicalVolume *, Doubles>(lv, Doubles(length, width)));
0192 }
0193 
0194 //________________________________________________________________________________________
0195 double DreamSD::curve_LY(const G4Step *aStep, int flag) {
0196   auto const stepPoint = aStep->GetPreStepPoint();
0197   auto const lv = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
0198   G4String nameVolume = lv->GetName();
0199 
0200   double weight = 1.;
0201   G4ThreeVector localPoint = setToLocal(stepPoint->GetPosition(), stepPoint->GetTouchable());
0202   double crlength = crystalLength(lv);
0203   double localz = localPoint.x();
0204   double dapd = 0.5 * crlength - flag * localz;  // Distance from closest APD
0205   if (dapd >= -0.1 || dapd <= crlength + 0.1) {
0206     if (dapd <= 100.)
0207       weight = 1.0 + slopeLY_ - dapd * 0.01 * slopeLY_;
0208   } else {
0209     edm::LogWarning("EcalSim") << "DreamSD: light coll curve : wrong distance "
0210                                << "to APD " << dapd << " crlength = " << crlength << " crystal name = " << nameVolume
0211                                << " z of localPoint = " << localz << " take weight = " << weight;
0212   }
0213 #ifdef EDM_ML_DEBUG
0214   edm::LogVerbatim("EcalSim") << "DreamSD, light coll curve : " << dapd << " crlength = " << crlength
0215                               << " crystal name = " << nameVolume << " z of localPoint = " << localz
0216                               << " take weight = " << weight;
0217 #endif
0218   return weight;
0219 }
0220 
0221 //________________________________________________________________________________________
0222 double DreamSD::crystalLength(G4LogicalVolume *lv) const {
0223   double length = -1.;
0224   DimensionMap::const_iterator ite = xtalLMap_.find(lv);
0225   if (ite != xtalLMap_.end())
0226     length = ite->second.first;
0227   return length;
0228 }
0229 
0230 //________________________________________________________________________________________
0231 double DreamSD::crystalWidth(G4LogicalVolume *lv) const {
0232   double width = -1.;
0233   DimensionMap::const_iterator ite = xtalLMap_.find(lv);
0234   if (ite != xtalLMap_.end())
0235     width = ite->second.second;
0236   return width;
0237 }
0238 
0239 //________________________________________________________________________________________
0240 // Calculate total cherenkov deposit
0241 // Inspired by Geant4's Cherenkov implementation
0242 double DreamSD::cherenkovDeposit_(const G4Step *aStep) {
0243   double cherenkovEnergy = 0;
0244   if (!materialPropertiesTable_)
0245     return cherenkovEnergy;
0246   G4Material *material = aStep->GetTrack()->GetMaterial();
0247 
0248   // Retrieve refractive index
0249   G4MaterialPropertyVector *Rindex = materialPropertiesTable_->GetProperty("RINDEX");
0250   if (Rindex == nullptr) {
0251     edm::LogWarning("EcalSim") << "Couldn't retrieve refractive index";
0252     return cherenkovEnergy;
0253   }
0254 
0255   // V.Ivanchenko - temporary close log output for 9.5
0256   // Material refraction properties
0257   int Rlength = Rindex->GetVectorLength() - 1;
0258   double Pmin = Rindex->Energy(0);
0259   double Pmax = Rindex->Energy(Rlength);
0260 #ifdef EDM_ML_DEBUG
0261   edm::LogVerbatim("EcalSim") << "Material properties: \n  Pmin = " << Pmin << "  Pmax = " << Pmax;
0262 #endif
0263   // Get particle properties
0264   const G4StepPoint *pPreStepPoint = aStep->GetPreStepPoint();
0265   const G4StepPoint *pPostStepPoint = aStep->GetPostStepPoint();
0266   const G4ThreeVector &x0 = pPreStepPoint->GetPosition();
0267   G4ThreeVector p0 = aStep->GetDeltaPosition().unit();
0268   const G4DynamicParticle *aParticle = aStep->GetTrack()->GetDynamicParticle();
0269   const double charge = aParticle->GetDefinition()->GetPDGCharge();
0270   // beta is averaged over step
0271   double beta = 0.5 * (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta());
0272   double BetaInverse = 1.0 / beta;
0273 
0274 #ifdef EDM_ML_DEBUG
0275   edm::LogVerbatim("EcalSim") << "Particle properties: \n  charge = " << charge << "  beta   = " << beta;
0276 #endif
0277 
0278   // Now get number of photons generated in this step
0279   double meanNumberOfPhotons = getAverageNumberOfPhotons_(charge, beta, material, Rindex);
0280   if (meanNumberOfPhotons <= 0.0) {  // Don't do anything
0281 #ifdef EDM_ML_DEBUG
0282     edm::LogVerbatim("EcalSim") << "Mean number of photons is zero: " << meanNumberOfPhotons << ", stopping here";
0283 #endif
0284     return cherenkovEnergy;
0285   }
0286 
0287   // number of photons is in unit of Geant4...
0288   meanNumberOfPhotons *= aStep->GetStepLength();
0289 
0290   // Now get a poisson distribution
0291   int numPhotons = static_cast<int>(G4Poisson(meanNumberOfPhotons));
0292   // edm::LogVerbatim("EcalSim") << "Number of photons = " << numPhotons;
0293   if (numPhotons <= 0) {
0294 #ifdef EDM_ML_DEBUG
0295     edm::LogVerbatim("EcalSim") << "Poission number of photons is zero: " << numPhotons << ", stopping here";
0296 #endif
0297     return cherenkovEnergy;
0298   }
0299 
0300   // Material refraction properties
0301   double dp = Pmax - Pmin;
0302   double maxCos = BetaInverse / (*Rindex)[Rlength];
0303   double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
0304 
0305   // Finally: get contribution of each photon
0306   for (int iPhoton = 0; iPhoton < numPhotons; ++iPhoton) {
0307     // Determine photon momentum
0308     double randomNumber;
0309     double sampledMomentum, sampledRI;
0310     double cosTheta, sin2Theta;
0311 
0312     // sample a momentum (not sure why this is needed!)
0313     do {
0314       randomNumber = G4UniformRand();
0315       sampledMomentum = Pmin + randomNumber * dp;
0316       sampledRI = Rindex->Value(sampledMomentum);
0317       cosTheta = BetaInverse / sampledRI;
0318 
0319       sin2Theta = (1.0 - cosTheta) * (1.0 + cosTheta);
0320       randomNumber = G4UniformRand();
0321 
0322     } while (randomNumber * maxSin2 > sin2Theta);
0323 
0324     // Generate random position of photon on cone surface
0325     // defined by Theta
0326     randomNumber = G4UniformRand();
0327 
0328     double phi = twopi * randomNumber;
0329     double sinPhi = sin(phi);
0330     double cosPhi = cos(phi);
0331 
0332     // Create photon momentum direction vector
0333     // The momentum direction is still w.r.t. the coordinate system where the
0334     // primary particle direction is aligned with the z axis
0335     double sinTheta = sqrt(sin2Theta);
0336     double px = sinTheta * cosPhi;
0337     double py = sinTheta * sinPhi;
0338     double pz = cosTheta;
0339     G4ThreeVector photonDirection(px, py, pz);
0340 
0341     // Rotate momentum direction back to global (crystal) reference system
0342     photonDirection.rotateUz(p0);
0343 
0344     // Create photon position and momentum
0345     randomNumber = G4UniformRand();
0346     G4ThreeVector photonPosition = x0 + randomNumber * aStep->GetDeltaPosition();
0347     G4ThreeVector photonMomentum = sampledMomentum * photonDirection;
0348 
0349     // Collect energy on APD
0350     cherenkovEnergy += getPhotonEnergyDeposit_(photonMomentum, photonPosition, aStep);
0351   }
0352   return cherenkovEnergy;
0353 }
0354 
0355 //________________________________________________________________________________________
0356 // Returns number of photons produced per GEANT-unit (millimeter) in the current
0357 // medium. From G4Cerenkov.cc
0358 double DreamSD::getAverageNumberOfPhotons_(const double charge,
0359                                            const double beta,
0360                                            const G4Material *aMaterial,
0361                                            const G4MaterialPropertyVector *Rindex) {
0362   const double rFact = 369.81 / (eV * cm);
0363 
0364   if (beta <= 0.0)
0365     return 0.0;
0366 
0367   double BetaInverse = 1. / beta;
0368 
0369   // Vectors used in computation of Cerenkov Angle Integral:
0370   //         - Refraction Indices for the current material
0371   //        - new G4PhysicsOrderedFreeVector allocated to hold CAI's
0372 
0373   // Min and Max photon momenta
0374   int Rlength = Rindex->GetVectorLength() - 1;
0375   double Pmin = Rindex->Energy(0);
0376   double Pmax = Rindex->Energy(Rlength);
0377 
0378   // Min and Max Refraction Indices
0379   double nMin = (*Rindex)[0];
0380   double nMax = (*Rindex)[Rlength];
0381 
0382   // Max Cerenkov Angle Integral
0383   double CAImax = chAngleIntegrals_.get()->GetMaxEnergy();
0384 
0385   double dp = 0., ge = 0., CAImin = 0.;
0386 
0387   // If n(Pmax) < 1/Beta -- no photons generated
0388   if (nMax < BetaInverse) {
0389   }
0390 
0391   // otherwise if n(Pmin) >= 1/Beta -- photons generated
0392   else if (nMin > BetaInverse) {
0393     dp = Pmax - Pmin;
0394     ge = CAImax;
0395   }
0396   // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
0397   // we need to find a P such that the value of n(P) == 1/Beta.
0398   // Interpolation is performed by the GetPhotonEnergy() and
0399   // GetProperty() methods of the G4MaterialPropertiesTable and
0400   // the GetValue() method of G4PhysicsVector.
0401   else {
0402     Pmin = Rindex->Value(BetaInverse);
0403     dp = Pmax - Pmin;
0404     // need boolean for current implementation of G4PhysicsVector
0405     // ==> being phased out
0406     double CAImin = chAngleIntegrals_->Value(Pmin);
0407     ge = CAImax - CAImin;
0408   }
0409 
0410   // Calculate number of photons
0411   double numPhotons = rFact * charge / eplus * charge / eplus * (dp - ge * BetaInverse * BetaInverse);
0412 
0413   edm::LogVerbatim("EcalSim") << "@SUB=getAverageNumberOfPhotons\nCAImin = " << CAImin << "\nCAImax = " << CAImax
0414                               << "\ndp = " << dp << ", ge = " << ge << "\nnumPhotons = " << numPhotons;
0415   return numPhotons;
0416 }
0417 
0418 //________________________________________________________________________________________
0419 // Set lead tungstate material properties on the fly.
0420 // Values from Ts42 detector construction
0421 bool DreamSD::setPbWO2MaterialProperties_(G4Material *aMaterial) {
0422   std::string pbWO2Name("E_PbWO4");
0423   std::string name = static_cast<std::string>(aMaterial->GetName());
0424   if (static_cast<std::string>(dd4hep::dd::noNamespace(name)) != pbWO2Name) {  // Wrong material!
0425     edm::LogWarning("EcalSim") << "This is not the right material: "
0426                                << "expecting " << pbWO2Name << ", got " << aMaterial->GetName();
0427     return false;
0428   }
0429 
0430   G4MaterialPropertiesTable *table = new G4MaterialPropertiesTable();
0431 
0432   // Refractive index as a function of photon momentum
0433   // FIXME: Should somehow put that in the configuration
0434   const int nEntries = 14;
0435   double PhotonEnergy[nEntries] = {1.7712 * eV,
0436                                    1.8368 * eV,
0437                                    1.90745 * eV,
0438                                    1.98375 * eV,
0439                                    2.0664 * eV,
0440                                    2.15625 * eV,
0441                                    2.25426 * eV,
0442                                    2.3616 * eV,
0443                                    2.47968 * eV,
0444                                    2.61019 * eV,
0445                                    2.75521 * eV,
0446                                    2.91728 * eV,
0447                                    3.09961 * eV,
0448                                    3.30625 * eV};
0449   double RefractiveIndex[nEntries] = {2.17728,
0450                                       2.18025,
0451                                       2.18357,
0452                                       2.18753,
0453                                       2.19285,
0454                                       2.19813,
0455                                       2.20441,
0456                                       2.21337,
0457                                       2.22328,
0458                                       2.23619,
0459                                       2.25203,
0460                                       2.27381,
0461                                       2.30282,
0462                                       2.34666};
0463 
0464   table->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex, nEntries);
0465   aMaterial->SetMaterialPropertiesTable(table);  // FIXME: could this leak? What does G4 do?
0466 
0467   // Calculate Cherenkov angle integrals:
0468   // This is an ad-hoc solution (we hold it in the class, not in the material)
0469   chAngleIntegrals_ = std::make_unique<G4PhysicsFreeVector>(nEntries);
0470 
0471   int index = 0;
0472   double currentRI = RefractiveIndex[index];
0473   double currentPM = PhotonEnergy[index];
0474   double currentCAI = 0.0;
0475   chAngleIntegrals_.get()->PutValue(0, currentPM, currentCAI);
0476   double prevPM = currentPM;
0477   double prevCAI = currentCAI;
0478   double prevRI = currentRI;
0479   while (++index < nEntries) {
0480     currentRI = RefractiveIndex[index];
0481     currentPM = PhotonEnergy[index];
0482     currentCAI = 0.5 * (1.0 / (prevRI * prevRI) + 1.0 / (currentRI * currentRI));
0483     currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
0484 
0485     chAngleIntegrals_.get()->PutValue(index, currentPM, currentCAI);
0486 
0487     prevPM = currentPM;
0488     prevCAI = currentCAI;
0489     prevRI = currentRI;
0490   }
0491 
0492 #ifdef EDM_ML_DEBUG
0493   edm::LogVerbatim("EcalSim") << "Material properties set for " << aMaterial->GetName();
0494 #endif
0495   return true;
0496 }
0497 
0498 //________________________________________________________________________________________
0499 // Calculate energy deposit of a photon on APD
0500 // - simple tracing to APD position (straight line);
0501 // - configurable reflection probability if not straight to APD;
0502 // - APD response function
0503 double DreamSD::getPhotonEnergyDeposit_(const G4ThreeVector &p, const G4ThreeVector &x, const G4Step *aStep) {
0504   double energy = 0;
0505 
0506   // Crystal dimensions
0507 
0508   // edm::LogVerbatim("EcalSim") << p << x;
0509 
0510   // 1. Check if this photon goes straight to the APD:
0511   //    - assume that APD is at x=xtalLength/2.0
0512   //    - extrapolate from x=x0 to x=xtalLength/2.0 using momentum in x-y
0513 
0514   G4StepPoint *stepPoint = aStep->GetPreStepPoint();
0515   G4LogicalVolume *lv = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
0516   G4String nameVolume = lv->GetName();
0517 
0518   double crlength = crystalLength(lv);
0519   double crwidth = crystalWidth(lv);
0520   double dapd = 0.5 * crlength - x.x();  // Distance from closest APD
0521   double y = p.y() / p.x() * dapd;
0522 
0523 #ifdef EDM_ML_DEBUG
0524   edm::LogVerbatim("EcalSim") << "Distance to APD: " << dapd << " - y at APD: " << y;
0525 #endif
0526   // Not straight: compute probability
0527   if (std::abs(y) > crwidth * 0.5) {
0528   }
0529 
0530   // 2. Retrieve efficiency for this wavelength (in nm, from MeV)
0531   double waveLength = p.mag() * 1.239e8;
0532 
0533   energy = p.mag() * PMTResponse::getEfficiency(waveLength);
0534 #ifdef EDM_ML_DEBUG
0535   edm::LogVerbatim("EcalSim") << "Wavelength: " << waveLength << " - Energy: " << energy;
0536 #endif
0537   return energy;
0538 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!