Par03/src/Par03EMShowerModel.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file Par03EMShowerModel.cc
0027 /// \brief Implementation of the Par03EMShowerModel class
0028
0029 #include "Par03EMShowerModel.hh"
0030
0031 #include "Par03EMShowerMessenger.hh"
0032
0033 #include "G4Electron.hh"
0034 #include "G4FastHit.hh"
0035 #include "G4FastSimHitMaker.hh"
0036 #include "G4Gamma.hh"
0037 #include "G4Positron.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4UnitsTable.hh"
0040 #include "Randomize.hh"
0041
0042 Par03EMShowerModel::Par03EMShowerModel(G4String aModelName, G4Region* aEnvelope)
0043   : G4VFastSimulationModel(aModelName, aEnvelope),
0044     fMessenger(new Par03EMShowerMessenger(this)),
0045     fHitMaker(new G4FastSimHitMaker)
0046 {}
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049
0050 Par03EMShowerModel::Par03EMShowerModel(G4String aModelName)
0051   : G4VFastSimulationModel(aModelName),
0052     fMessenger(new Par03EMShowerMessenger(this)),
0053     fHitMaker(new G4FastSimHitMaker)
0054 {}
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 Par03EMShowerModel::~Par03EMShowerModel() = default;
0059
0060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0061
0062 G4bool Par03EMShowerModel::IsApplicable(const G4ParticleDefinition& aParticleType)
0063 {
0064   return &aParticleType == G4Electron::ElectronDefinition()
0065          || &aParticleType == G4Positron::PositronDefinition()
0066          || &aParticleType == G4Gamma::GammaDefinition();
0067 }
0068
0069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0070
0071 G4bool Par03EMShowerModel::ModelTrigger(const G4FastTrack& aFastTrack)
0072 {
0073   // Check energy
0074   if (aFastTrack.GetPrimaryTrack()->GetKineticEnergy() < 1 * GeV) {
0075     return false;
0076   }
0077   // Check length of detector
0078   // Calculate depth of the detector along shower axis to verify if shower
0079   // will fit inside. Required max shower depth is defined by fLongMaxDepth, and
0080   // can be changed with UI command `/Par03/fastSim/longitudinalProfile/maxDepth
0081   G4double X0 = aFastTrack.GetPrimaryTrack()->GetMaterial()->GetRadlen();
0082   auto particleDirection = aFastTrack.GetPrimaryTrackLocalDirection();
0083   auto particlePosition = aFastTrack.GetPrimaryTrackLocalPosition();
0084   G4double detectorDepthInMM =
0085     aFastTrack.GetEnvelopeSolid()->DistanceToOut(particlePosition, particleDirection);
0086   G4double detectorDepthInX0 = detectorDepthInMM / X0;
0087   // check if detector depth is sufficient to create showers
0088   if (detectorDepthInX0 < fLongMaxDepth) {
0089     return false;
0090   }
0091   return true;
0092 }
0093
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095
0096 void Par03EMShowerModel::DoIt(const G4FastTrack& aFastTrack, G4FastStep& aFastStep)
0097 {
0098   // Remove particle from further processing by G4
0099   aFastStep.KillPrimaryTrack();
0100   aFastStep.ProposePrimaryTrackPathLength(0.0);
0101   G4double energy = aFastTrack.GetPrimaryTrack()->GetKineticEnergy();
0102   // No need to create any deposit, it will be handled by this model (and
0103   // G4FastSimHitMaker that will call the sensitive detector)
0104   aFastStep.ProposeTotalEnergyDeposited(0);
0105   auto particlePosition = aFastTrack.GetPrimaryTrackLocalPosition();
0106   auto particleDirection = aFastTrack.GetPrimaryTrackLocalDirection();
0107
0108   // Calculate how to create energy deposits
0109   // Following PDG 33.5 chapter
0110   // material calculation assumes homogeneous detector (true for Par03 example)
0111   auto material = aFastTrack.GetPrimaryTrack()->GetMaterial();
0112   G4double materialX0 = material->GetRadlen();
0113   G4double materialZ = material->GetZ();
0114   // EC estimation follows PDG fit to solids in Fig. 33.14 (rms 2.2%)
0115   G4double materialEc = 610 * MeV / (materialZ + 1.24);
0116   // RM estimation follows PDG Eq. (33.37) (rms 2.2%)
0117   G4double materialRM = 21.2052 * MeV * materialX0 / materialEc;
0118   G4double particleY = energy / materialEc;
0119   // Estimate shower maximum and alpha parameter of Gamma distribution
0120   // that describes the longitudinal profile (PDG Eq. (33.35))
0121   // unless alpha is specified by UI command
0122   if (fAlpha < 0) {
0123     // from PDG Eq. (33.36)
0124     G4double particleTmax = std::log(particleY);
0125     if (aFastTrack.GetPrimaryTrack()->GetParticleDefinition() == G4Gamma::GammaDefinition()) {
0126       particleTmax += 0.5;
0127     }
0128     else {
0129       particleTmax -= 0.5;
0130     }
0131     fAlpha = particleTmax * fBeta + 1;
0132   }
0133   // Unless sigma of Gaussian distribution describing the transverse profile
0134   // is specified by UI command, use value calculated from Moliere Radius
0135   if (fSigma < 0) {
0136     // 90% of shower is contained within 1 * R_M
0137     // 1.645 * std dev of Gaussian contains 90%
0138     fSigma = materialRM / 1.645;
0139   }
0140
0141   // Calculate rotation matrix along the particle momentum direction
0142   // It will rotate the shower axes to match the incoming particle direction
0143   G4RotationMatrix rotMatrix = G4RotationMatrix();
0144   double particleTheta = particleDirection.theta();
0145   double particlePhi = particleDirection.phi();
0146   double epsilon = 1e-3;
0147   rotMatrix.rotateY(particleTheta);
0148   // do not use (random) phi if x==y==0
0149   if (!(std::fabs(particleDirection.x()) < epsilon && std::fabs(particleDirection.y()) < epsilon))
0150     rotMatrix.rotateZ(particlePhi);
0151
0152   // Create hits
0153   // First use rejecton sampling to sample from Gamma distribution
0154   // then get random numbers from uniform distribution for azimuthal angle, and
0155   // from Gaussian for radius
0156   G4ThreeVector position;
0157   G4double gammaMax = Gamma((fAlpha - 1) / fBeta, fAlpha, fBeta);
0158   G4int generatedHits = 0;
0159   while (generatedHits < fNbOfHits) {
0160     G4double random1 = G4UniformRand() * fLongMaxDepth;
0161     G4double random2 = G4UniformRand() * gammaMax;
0162     if (Gamma(random1, fAlpha, fBeta) >= random2) {
0163       // Generate corresponding rho (phi) from Gaussian (flat) distribution
0164       G4double phiPosition = G4UniformRand() * 2 * CLHEP::pi;
0165       G4double rhoPosition = G4RandGauss::shoot(0, fSigma);
0166       position = particlePosition
0167                  + rotMatrix
0168                      * G4ThreeVector(rhoPosition * std::sin(phiPosition),
0169                                      rhoPosition * std::cos(phiPosition), random1 * materialX0);
0170       // Create energy deposit in the detector
0171       // This will call appropriate sensitive detector class
0172       fHitMaker->make(G4FastHit(position, energy / fNbOfHits), aFastTrack);
0173       generatedHits++;
0174     }
0175   }
0176 }
0177
0178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0179
0180 void Par03EMShowerModel::Print() const
0181 {
0182   G4cout << "Par03EMShowerModel: " << G4endl;
0183   G4cout << "Gaussian distribution (transverse plane): \tmu = 0, sigma = "
0184          << G4BestUnit(fSigma, "Length") << G4endl;
0185   if (fSigma < 0)
0186     G4cout << "Negative sigma value means that it will be recalculated "
0187               "from the value of the Moliere radius of the detector material, "
0188               "taking into account that 90% of the area below the Gaussian "
0189               "distribution (from mu - 1.645 sigma to mu + 1.645 sigma) "
0190               "corresponds to area within 1 Moliere radius."
0191            << G4endl;
0192   G4cout << "Gamma distribution (along shower axis): \talpha = " << fAlpha << ", beta = " << fBeta
0193          << ", max depth = " << fLongMaxDepth << " X0" << G4endl;
0194   if (fAlpha < 0)
0195     G4cout << "Negative alpha value means that it will be recalculated "
0196               "from the critical energy of the detector material, particle "
0197               "type, and beta parameter.\n alpha = beta * T_max, where T_max = "
0198               "ln(E/E_C) + C\n where E is particle energy, E_C is critical "
0199               "energy in the material, and constant C = -0.5 for electrons and "
0200               "0.5 for photons (Eq. (33.36) from PDG)."
0201            << G4endl;
0202   G4cout << "Number of created energy deposits: " << fNbOfHits << G4endl;
0203 }