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
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020
0021
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
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
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
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
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
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
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;
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
0241
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
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
0256
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
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
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
0279 double meanNumberOfPhotons = getAverageNumberOfPhotons_(charge, beta, material, Rindex);
0280 if (meanNumberOfPhotons <= 0.0) {
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
0288 meanNumberOfPhotons *= aStep->GetStepLength();
0289
0290
0291 int numPhotons = static_cast<int>(G4Poisson(meanNumberOfPhotons));
0292
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
0301 double dp = Pmax - Pmin;
0302 double maxCos = BetaInverse / (*Rindex)[Rlength];
0303 double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
0304
0305
0306 for (int iPhoton = 0; iPhoton < numPhotons; ++iPhoton) {
0307
0308 double randomNumber;
0309 double sampledMomentum, sampledRI;
0310 double cosTheta, sin2Theta;
0311
0312
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
0325
0326 randomNumber = G4UniformRand();
0327
0328 double phi = twopi * randomNumber;
0329 double sinPhi = sin(phi);
0330 double cosPhi = cos(phi);
0331
0332
0333
0334
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
0342 photonDirection.rotateUz(p0);
0343
0344
0345 randomNumber = G4UniformRand();
0346 G4ThreeVector photonPosition = x0 + randomNumber * aStep->GetDeltaPosition();
0347 G4ThreeVector photonMomentum = sampledMomentum * photonDirection;
0348
0349
0350 cherenkovEnergy += getPhotonEnergyDeposit_(photonMomentum, photonPosition, aStep);
0351 }
0352 return cherenkovEnergy;
0353 }
0354
0355
0356
0357
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
0370
0371
0372
0373
0374 int Rlength = Rindex->GetVectorLength() - 1;
0375 double Pmin = Rindex->Energy(0);
0376 double Pmax = Rindex->Energy(Rlength);
0377
0378
0379 double nMin = (*Rindex)[0];
0380 double nMax = (*Rindex)[Rlength];
0381
0382
0383 double CAImax = chAngleIntegrals_.get()->GetMaxEnergy();
0384
0385 double dp = 0., ge = 0., CAImin = 0.;
0386
0387
0388 if (nMax < BetaInverse) {
0389 }
0390
0391
0392 else if (nMin > BetaInverse) {
0393 dp = Pmax - Pmin;
0394 ge = CAImax;
0395 }
0396
0397
0398
0399
0400
0401 else {
0402 Pmin = Rindex->Value(BetaInverse);
0403 dp = Pmax - Pmin;
0404
0405
0406 double CAImin = chAngleIntegrals_->Value(Pmin);
0407 ge = CAImax - CAImin;
0408 }
0409
0410
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
0420
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) {
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
0433
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);
0466
0467
0468
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
0500
0501
0502
0503 double DreamSD::getPhotonEnergyDeposit_(const G4ThreeVector &p, const G4ThreeVector &x, const G4Step *aStep) {
0504 double energy = 0;
0505
0506
0507
0508
0509
0510
0511
0512
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();
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
0527 if (std::abs(y) > crwidth * 0.5) {
0528 }
0529
0530
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 }