Application/src/LowEnergyFastSimModel.cc

0001 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "SimG4Core/Application/interface/LowEnergyFastSimModel.h"
0005 #include "SimG4Core/Application/interface/TrackingAction.h"
0006 #include "SimG4Core/Geometry/interface/DD4hep2DDDName.h"
0007 #include "SimG4Core/Notification/interface/TrackInformation.h"
0008
0009 #include "G4VFastSimulationModel.hh"
0010 #include "G4EventManager.hh"
0011 #include "G4Electron.hh"
0012 #include "GFlashHitMaker.hh"
0013 #include "G4Region.hh"
0014 #include "G4Material.hh"
0015 #include "G4Positron.hh"
0016 #include "G4ParticleDefinition.hh"
0017 #include "G4PhysicalConstants.hh"
0018
0019 constexpr G4double twomass = 2 * CLHEP::electron_mass_c2;
0020 constexpr G4double scaleFactor = 1.025;
0021
0022 LowEnergyFastSimModel::LowEnergyFastSimModel(const G4String& name, G4Region* region, const edm::ParameterSet& parSet)
0023     : G4VFastSimulationModel(name, region),
0024       fRegion(region),
0025       fTrackingAction(nullptr),
0026       fCheck(false),
0027       fTailPos(0., 0., 0.) {
0028   fEmax = parSet.getParameter<double>("LowEnergyGflashEcalEmax") * CLHEP::GeV;
0029   fPositron = G4Positron::Positron();
0030   fMaterial = nullptr;
0031   auto table = G4Material::GetMaterialTable();
0032   for (auto const& mat : *table) {
0033     const G4String& nam = mat->GetName();
0034     std::size_t n = nam.size();
0035     if (n > 4) {
0036       const G4String& sn = nam.substr(n - 5, 5);
0037       if (sn == "PbWO4") {
0038         fMaterial = mat;
0039         break;
0040       }
0041     }
0042   }
0043   const G4String& nm = (nullptr == fMaterial) ? "not found!" : fMaterial->GetName();
0044   edm::LogVerbatim("LowEnergyFastSimModel") << "LowEGFlash material: <" << nm << ">";
0045 }
0046
0047 G4bool LowEnergyFastSimModel::IsApplicable(const G4ParticleDefinition& particle) {
0048   return (11 == std::abs(particle.GetPDGEncoding()));
0049 }
0050
0051 G4bool LowEnergyFastSimModel::ModelTrigger(const G4FastTrack& fastTrack) {
0052   const G4Track* track = fastTrack.GetPrimaryTrack();
0053   G4double energy = track->GetKineticEnergy();
0054   if (fMaterial != track->GetMaterial() || energy >= fEmax)
0055     return false;
0056
0057   /*
0058   edm::LogVerbatim("LowEnergyFastSimModel") << track->GetDefinition()->GetParticleName()
0059                         << " Ekin(MeV)=" << energy << " material: <"
0060                                             << track->GetMaterial()->GetName() << ">";
0061   */
0062   if (fCheck) {
0063     if (nullptr == fTrackingAction) {
0064       fTrackingAction = static_cast<const TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
0065     }
0066     const TrackInformation* ptrti = static_cast<TrackInformation*>(track->GetUserInformation());
0067     int pdg = ptrti->genParticlePID();
0068     if (std::abs(pdg) == 11 || pdg == 22)
0069       return false;
0070   }
0071   return true;
0072 }
0073
0074 void LowEnergyFastSimModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) {
0075   auto track = fastTrack.GetPrimaryTrack();
0076   G4double energy = track->GetKineticEnergy() * scaleFactor;
0077
0078   const G4ThreeVector& pos = track->GetPosition();
0079
0080   G4double inPointEnergy = fParam.GetInPointEnergyFraction(energy) * energy;
0081
0082   // take into account positron annihilation (not included in in-point)
0083   if (fPositron == track->GetDefinition())
0084     energy += twomass;
0085
0086   const G4ThreeVector& momDir = track->GetMomentumDirection();
0087
0088   // Russian roulette
0089   double wt2 = track->GetWeight();
0090   if (wt2 <= 0.0) {
0091     wt2 = 1.0;
0092   }
0093   // tail energy deposition
0094   G4double etail = energy - inPointEnergy;
0095   G4int nspots = G4lrint(etail);
0096   if (0 == nspots) {
0097     inPointEnergy = energy;
0098     etail = 0.0;
0099   }
0100
0101   // in point energy deposition
0102   GFlashEnergySpot spot;
0103   spot.SetEnergy(inPointEnergy * wt2);
0104   spot.SetPosition(pos);
0105   fHitMaker.make(&spot, &fastTrack);
0106
0107   G4double zz = 0.0;
0108   if (0 < nspots) {
0109     etail *= wt2 / (G4double)nspots;
0110     /*
0111   edm::LogVerbatim("LowEnergyFastSimModel") << track->GetDefinition()->GetParticleName()
0112                         << " Ekin(MeV)=" << energy << " material: <"
0113                                             << track->GetMaterial()->GetName()
0114                         << "> Elocal=" << inPointEnergy
0115                         << " Etail=" << tailEnergy
0116                         << " Nspots=" << nspots;
0117   */
0118     for (G4int i = 0; i < nspots; ++i) {
0119       const G4double r = fParam.GetRadius(energy);
0120       const G4double z = fParam.GetZ();
0121       zz += z;
0122
0123       const G4double phi = CLHEP::twopi * G4UniformRand();
0124       fTailPos.set(r * std::cos(phi), r * std::sin(phi), z);
0125       fTailPos.rotateUz(momDir);
0126       fTailPos += pos;
0127
0128       spot.SetEnergy(etail);
0129       spot.SetPosition(fTailPos);
0130       fHitMaker.make(&spot, &fastTrack);
0131     }
0132     zz /= (G4double)nspots;
0133   }
0134   fastStep.KillPrimaryTrack();
0135   fastStep.ProposePrimaryTrackPathLength(zz);
0136 }