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(true),
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     G4String nam = (G4String)(DD4hep2DDDName::noNameSpace(mat->GetName()));
0034     size_t n = nam.size();
0035     if (n > 4) {
0036       G4String sn = nam.substr(n - 5, 5);
0037       if (sn == "PbWO4") {
0038         fMaterial = mat;
0039         break;
0040       }
0041     }
0042   }
0043   G4String nm = (nullptr == fMaterial) ? "not found!" : (G4String)(DD4hep2DDDName::noNameSpace(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   fastStep.KillPrimaryTrack();
0076   fastStep.ProposePrimaryTrackPathLength(0.0);
0077   auto track = fastTrack.GetPrimaryTrack();
0078   G4double energy = track->GetKineticEnergy() * scaleFactor;
0079
0080   const G4ThreeVector& pos = track->GetPosition();
0081
0082   G4double inPointEnergy = fParam.GetInPointEnergyFraction(energy) * energy;
0083
0084   // take into account positron annihilation (not included in in-point)
0085   if (fPositron == track->GetDefinition())
0086     energy += twomass;
0087
0088   const G4ThreeVector& momDir = track->GetMomentumDirection();
0089
0090   // in point energy deposition
0091   GFlashEnergySpot spot;
0092   spot.SetEnergy(inPointEnergy);
0093   spot.SetPosition(pos);
0094   fHitMaker.make(&spot, &fastTrack);
0095
0096   // Russian roulette
0097   double wt2 = track->GetWeight();
0098   if (wt2 <= 0.0) {
0099     wt2 = 1.0;
0100   }
0101
0102   // tail energy deposition
0103   const G4double etail = energy - inPointEnergy;
0104   const G4int nspots = etail;
0105   const G4double tailEnergy = etail * wt2 / (nspots + 1);
0106   /*
0107   edm::LogVerbatim("LowEnergyFastSimModel") << track->GetDefinition()->GetParticleName()
0108                         << " Ekin(MeV)=" << energy << " material: <"
0109                                             << track->GetMaterial()->GetName()
0110                         << "> Elocal=" << inPointEnergy
0111                         << " Etail=" << tailEnergy
0112                         << " Nspots=" << nspots+1;
0113   */
0114   for (G4int i = 0; i <= nspots; ++i) {
0115     const G4double r = fParam.GetRadius(energy);
0116     const G4double z = fParam.GetZ();
0117
0118     const G4double phi = CLHEP::twopi * G4UniformRand();
0119     fTailPos.set(r * std::cos(phi), r * std::sin(phi), z);
0120     fTailPos.rotateUz(momDir);
0121     fTailPos += pos;
0122
0123     spot.SetEnergy(tailEnergy);
0124     spot.SetPosition(fTailPos);
0125     fHitMaker.make(&spot, &fastTrack);
0126   }
0127 }