Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:18 

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