File indexing completed on 2022-12-12 23:53:08
0001
0002
0003
0004
0005
0006 #include "SimG4CMS/Calo/interface/HFShower.h"
0007 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
0008
0009 #include "G4NavigationHistory.hh"
0010 #include "G4Step.hh"
0011 #include "G4Track.hh"
0012 #include "G4VSolid.hh"
0013 #include "Randomize.hh"
0014 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0015 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0016
0017
0018
0019 #include <iostream>
0020
0021 HFShower::HFShower(const std::string &name,
0022 const HcalDDDSimConstants *hcons,
0023 const HcalSimulationParameters *hps,
0024 edm::ParameterSet const &p,
0025 int chk)
0026 : cherenkov_(p.getParameter<edm::ParameterSet>("HFShower")), fibre_(hcons, hps, p), chkFibre_(chk) {
0027 edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
0028 applyFidCut_ = m_HF.getParameter<bool>("ApplyFiducialCut");
0029 edm::ParameterSet m_HF2 = m_HF.getParameter<edm::ParameterSet>("HFShowerBlock");
0030 equalizeTimeShift_ = m_HF2.getParameter<bool>("EqualizeTimeShift");
0031 probMax_ = m_HF2.getParameter<double>("ProbMax");
0032
0033 edm::LogVerbatim("HFShower") << "HFShower:: Maximum probability cut off " << probMax_ << " Check flag " << chkFibre_
0034 << " EqualizeTimeShift " << equalizeTimeShift_;
0035
0036
0037 gpar_ = hcons->getGparHF();
0038 }
0039
0040 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, double weight) {
0041 std::vector<HFShower::Hit> hits;
0042 int nHit = 0;
0043
0044 double edep = weight * (aStep->GetTotalEnergyDeposit());
0045 #ifdef EDM_ML_DEBUG
0046 edm::LogVerbatim("HFShower") << "HFShower::getHits: energy " << aStep->GetTotalEnergyDeposit() << " weight " << weight
0047 << " edep " << edep;
0048 #endif
0049 double stepl = 0.;
0050
0051 if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
0052 stepl = aStep->GetStepLength();
0053 if ((edep == 0.) || (stepl == 0.)) {
0054 #ifdef EDM_ML_DEBUG
0055 edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0056 #endif
0057 return hits;
0058 }
0059 const G4Track *aTrack = aStep->GetTrack();
0060 const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
0061
0062 HFShower::Hit hit;
0063 double energy = aParticle->GetTotalEnergy();
0064 double momentum = aParticle->GetTotalMomentum();
0065 double pBeta = momentum / energy;
0066 double dose = 0.;
0067 int npeDose = 0;
0068
0069 const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
0070 const G4ParticleDefinition *particleDef = aTrack->GetDefinition();
0071
0072 const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0073 const G4ThreeVector &globalPos = preStepPoint->GetPosition();
0074 G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
0075
0076 double zv = std::abs(globalPos.z()) - gpar_[4];
0077 G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
0078 G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
0079
0080 int depth = 1;
0081 int npmt = 0;
0082 bool ok = true;
0083 if (zv < 0. || zv > gpar_[1]) {
0084 ok = false;
0085 }
0086 if (ok && applyFidCut_) {
0087 npmt = HFFibreFiducial::PMTNumber(globalPos);
0088 if (npmt <= 0) {
0089 ok = false;
0090 } else if (npmt > 24) {
0091 if (zv > gpar_[0]) {
0092 depth = 2;
0093 } else {
0094 ok = false;
0095 }
0096 }
0097 #ifdef EDM_ML_DEBUG
0098 edm::LogVerbatim("HFShower") << "HFShower: npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":" << gpar_[4]
0099 << ":" << zv << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
0100 #endif
0101 } else {
0102 depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10;
0103 }
0104 G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
0105
0106 double u = localMom.x();
0107 double v = localMom.y();
0108 double w = localMom.z();
0109 double zCoor = localPos.z();
0110 double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
0111 double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
0112 double time =
0113 (equalizeTimeShift_) ? (fibre_.tShift(localPos, depth, -1)) : (fibre_.tShift(localPos, depth, chkFibre_));
0114
0115 #ifdef EDM_ML_DEBUG
0116 edm::LogVerbatim("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
0117 << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
0118 << "\n Direction " << momentumDir << " Local " << localMom;
0119 #endif
0120
0121 int npe = 1;
0122 std::vector<double> wavelength;
0123 std::vector<double> momz;
0124 if (!applyFidCut_) {
0125 if (ok)
0126 npe = cherenkov_.computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
0127 wavelength = cherenkov_.getWL();
0128 momz = cherenkov_.getMom();
0129 }
0130 if (ok && npe > 0) {
0131 for (int i = 0; i < npe; ++i) {
0132 double p = 1.;
0133 if (!applyFidCut_)
0134 p = fibre_.attLength(wavelength[i]);
0135 double r1 = G4UniformRand();
0136 double r2 = G4UniformRand();
0137 #ifdef EDM_ML_DEBUG
0138 edm::LogVerbatim("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1 << ":" << exp(-p * zFibre)
0139 << " r2 " << r2 << ":" << probMax_
0140 << " Survive: " << (r1 <= exp(-p * zFibre) && r2 <= probMax_);
0141 #endif
0142 if (applyFidCut_ || chkFibre_ < 0 || (r1 <= exp(-p * zFibre) && r2 <= probMax_)) {
0143 hit.depth = depth;
0144 hit.time = tSlice + time;
0145
0146 if (!applyFidCut_) {
0147 hit.wavelength = wavelength[i];
0148 hit.momentum = momz[i];
0149 } else {
0150 hit.wavelength = 300.;
0151 hit.momentum = 1.;
0152 }
0153 hit.position = globalPos;
0154 hits.push_back(hit);
0155 nHit++;
0156 }
0157 }
0158 }
0159
0160 #ifdef EDM_ML_DEBUG
0161 edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0162 for (int i = 0; i < nHit; ++i)
0163 edm::LogVerbatim("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time "
0164 << hits[i].time << " Momentum " << hits[i].momentum << " Position "
0165 << hits[i].position;
0166 #endif
0167 return hits;
0168 }
0169
0170 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, bool forLibraryProducer, double zoffset) {
0171 std::vector<HFShower::Hit> hits;
0172 int nHit = 0;
0173
0174 double edep = aStep->GetTotalEnergyDeposit();
0175 double stepl = 0.;
0176
0177 if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
0178 stepl = aStep->GetStepLength();
0179 if ((edep == 0.) || (stepl == 0.)) {
0180 #ifdef EDM_ML_DEBUG
0181 edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0182 #endif
0183 return hits;
0184 }
0185 const G4Track *aTrack = aStep->GetTrack();
0186 const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
0187
0188 HFShower::Hit hit;
0189 double energy = aParticle->GetTotalEnergy();
0190 double momentum = aParticle->GetTotalMomentum();
0191 double pBeta = momentum / energy;
0192 double dose = 0.;
0193 int npeDose = 0;
0194
0195 const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
0196 G4ParticleDefinition *particleDef = aTrack->GetDefinition();
0197
0198 G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0199 const G4ThreeVector &globalPos = preStepPoint->GetPosition();
0200 G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
0201 double zb = std::abs(globalPos.z() - zoffset);
0202 double zv = zb - .5 * gpar_[1];
0203 G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
0204 G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
0205 bool ok = true;
0206 if (zb < 0. || zb > gpar_[1]) {
0207 ok = false;
0208 }
0209 int depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10;
0210 G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
0211
0212 double u = localMom.x();
0213 double v = localMom.y();
0214 double w = localMom.z();
0215 double zCoor = localPos.z();
0216 double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
0217 double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
0218
0219 #ifdef EDM_ML_DEBUG
0220 double time = (equalizeTimeShift_) ? 0 : fibre_.tShift(localPos, depth, chkFibre_);
0221 edm::LogVerbatim("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
0222 << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
0223 << "\n Direction " << momentumDir << " Local " << localMom;
0224 #endif
0225
0226 int npe = 0;
0227 std::vector<double> wavelength;
0228 std::vector<double> momz;
0229 if (ok)
0230 npe = cherenkov_.computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
0231 wavelength = cherenkov_.getWL();
0232 momz = cherenkov_.getMom();
0233 if (ok && npe > 0) {
0234 for (int i = 0; i < npe; ++i) {
0235 hit.depth = depth;
0236
0237 hit.time = tSlice;
0238 hit.wavelength = wavelength[i];
0239 hit.momentum = momz[i];
0240 hit.position = globalPos;
0241 hits.push_back(hit);
0242 nHit++;
0243 }
0244 }
0245
0246 #ifdef EDM_ML_DEBUG
0247 edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0248 for (int i = 0; i < nHit; ++i)
0249 edm::LogVerbatim("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time "
0250 << hits[i].time << " Momentum " << hits[i].momentum << " Position "
0251 << hits[i].position;
0252 #endif
0253 return hits;
0254 }
0255
0256 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, bool forLibrary) {
0257 std::vector<HFShower::Hit> hits;
0258 int nHit = 0;
0259
0260 double edep = aStep->GetTotalEnergyDeposit();
0261 double stepl = 0.;
0262
0263 if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
0264 stepl = aStep->GetStepLength();
0265 if ((edep == 0.) || (stepl == 0.)) {
0266 #ifdef EDM_ML_DEBUG
0267 edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0268 #endif
0269 return hits;
0270 }
0271
0272 const G4Track *aTrack = aStep->GetTrack();
0273 const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
0274
0275 HFShower::Hit hit;
0276 double energy = aParticle->GetTotalEnergy();
0277 double momentum = aParticle->GetTotalMomentum();
0278 double pBeta = momentum / energy;
0279 double dose = 0.;
0280 int npeDose = 0;
0281
0282 const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
0283 G4ParticleDefinition *particleDef = aTrack->GetDefinition();
0284
0285 const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0286 const G4ThreeVector &globalPos = preStepPoint->GetPosition();
0287 G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
0288 double zv = std::abs(globalPos.z()) - gpar_[4] - 0.5 * gpar_[1];
0289 G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
0290 G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
0291
0292 int depth = 1;
0293 int npmt = 0;
0294 bool ok = true;
0295 if (zv < 0 || zv > gpar_[1]) {
0296 ok = false;
0297 }
0298 if (ok && applyFidCut_) {
0299 npmt = HFFibreFiducial::PMTNumber(globalPos);
0300 if (npmt <= 0) {
0301 ok = false;
0302 } else if (npmt > 24) {
0303 if (zv > gpar_[0]) {
0304 depth = 2;
0305 } else {
0306 ok = false;
0307 }
0308 }
0309 #ifdef EDM_ML_DEBUG
0310 edm::LogVerbatim("HFShower") << "HFShower:getHits(SL): npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":"
0311 << gpar_[4] << ":" << zv << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
0312 #endif
0313 } else {
0314 depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10;
0315 }
0316 G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
0317
0318 double u = localMom.x();
0319 double v = localMom.y();
0320 double w = localMom.z();
0321 double zCoor = localPos.z();
0322 double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
0323 double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
0324 double time = (equalizeTimeShift_) ? 0 : fibre_.tShift(localPos, depth, chkFibre_);
0325
0326 #ifdef EDM_ML_DEBUG
0327 edm::LogVerbatim("HFShower") << "HFShower::getHits(SL): in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
0328 << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
0329 << "\n Direction " << momentumDir << " Local " << localMom;
0330 #endif
0331
0332 int npe = 0;
0333 if (ok)
0334 npe = cherenkov_.computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
0335 std::vector<double> wavelength = cherenkov_.getWL();
0336 std::vector<double> momz = cherenkov_.getMom();
0337
0338 for (int i = 0; i < npe; ++i) {
0339 hit.depth = depth;
0340 hit.time = tSlice + time;
0341 hit.wavelength = wavelength[i];
0342 hit.momentum = momz[i];
0343 hit.position = globalPos;
0344 hits.push_back(hit);
0345 nHit++;
0346 }
0347
0348 return hits;
0349 }