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