| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 #include "Par03EMShowerModel.hh"
27
28 #include "Par03EMShowerMessenger.hh"
29
30 #include "G4Electron.hh"
31 #include "G4FastHit.hh"
32 #include "G4FastSimHitMaker.hh"
33 #include "G4Gamma.hh"
34 #include "G4Positron.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4UnitsTable.hh"
37 #include "Randomize.hh"
38
39 Par03EMShowerModel::Par03EMShowerModel(G4String aModelName, G4Region* aEnvelope)
40 : G4VFastSimulationModel(aModelName, aEnvelope),
41 fMessenger(new Par03EMShowerMessenger(this)),
42 fHitMaker(new G4FastSimHitMaker)
43 {}
44
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46
47 Par03EMShowerModel::Par03EMShowerModel(G4String aModelName)
48 : G4VFastSimulationModel(aModelName),
49 fMessenger(new Par03EMShowerMessenger(this)),
50 fHitMaker(new G4FastSimHitMaker)
51 {}
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54
55 Par03EMShowerModel::~Par03EMShowerModel() = default;
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59 G4bool Par03EMShowerModel::IsApplicable(const G4ParticleDefinition& aParticleType)
60 {
61 return &aParticleType == G4Electron::ElectronDefinition()
62 || &aParticleType == G4Positron::PositronDefinition()
63 || &aParticleType == G4Gamma::GammaDefinition();
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67
68 G4bool Par03EMShowerModel::ModelTrigger(const G4FastTrack& aFastTrack)
69 {
70 // Check energy
71 if (aFastTrack.GetPrimaryTrack()->GetKineticEnergy() < 1 * GeV) {
72 return false;
73 }
74 // Check length of detector
75 // Calculate depth of the detector along shower axis to verify if shower
76 // will fit inside. Required max shower depth is defined by fLongMaxDepth, and
77 // can be changed with UI command `/Par03/fastSim/longitudinalProfile/maxDepth
78 G4double X0 = aFastTrack.GetPrimaryTrack()->GetMaterial()->GetRadlen();
79 auto particleDirection = aFastTrack.GetPrimaryTrackLocalDirection();
80 auto particlePosition = aFastTrack.GetPrimaryTrackLocalPosition();
81 G4double detectorDepthInMM =
82 aFastTrack.GetEnvelopeSolid()->DistanceToOut(particlePosition, particleDirection);
83 G4double detectorDepthInX0 = detectorDepthInMM / X0;
84 // check if detector depth is sufficient to create showers
85 if (detectorDepthInX0 < fLongMaxDepth) {
86 return false;
87 }
88 return true;
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92
93 void Par03EMShowerModel::DoIt(const G4FastTrack& aFastTrack, G4FastStep& aFastStep)
94 {
95 // Remove particle from further processing by G4
96 aFastStep.KillPrimaryTrack();
97 aFastStep.ProposePrimaryTrackPathLength(0.0);
98 G4double energy = aFastTrack.GetPrimaryTrack()->GetKineticEnergy();
99 // No need to create any deposit, it will be handled by this model (and
100 // G4FastSimHitMaker that will call the sensitive detector)
101 aFastStep.ProposeTotalEnergyDeposited(0);
102 auto particlePosition = aFastTrack.GetPrimaryTrackLocalPosition();
103 auto particleDirection = aFastTrack.GetPrimaryTrackLocalDirection();
104
105 // Calculate how to create energy deposits
106 // Following PDG 33.5 chapter
107 // material calculation assumes homogeneous detector (true for Par03 example)
108 auto material = aFastTrack.GetPrimaryTrack()->GetMaterial();
109 G4double materialX0 = material->GetRadlen();
110 G4double materialZ = material->GetZ();
111 // EC estimation follows PDG fit to solids in Fig. 33.14 (rms 2.2%)
112 G4double materialEc = 610 * MeV / (materialZ + 1.24);
113 // RM estimation follows PDG Eq. (33.37) (rms 2.2%)
114 G4double materialRM = 21.2052 * MeV * materialX0 / materialEc;
115 G4double particleY = energy / materialEc;
116 // Estimate shower maximum and alpha parameter of Gamma distribution
117 // that describes the longitudinal profile (PDG Eq. (33.35))
118 // unless alpha is specified by UI command
119 if (fAlpha < 0) {
120 // from PDG Eq. (33.36)
121 G4double particleTmax = std::log(particleY);
122 if (aFastTrack.GetPrimaryTrack()->GetParticleDefinition() == G4Gamma::GammaDefinition()) {
123 particleTmax += 0.5;
124 }
125 else {
126 particleTmax -= 0.5;
127 }
128 fAlpha = particleTmax * fBeta + 1;
129 }
130 // Unless sigma of Gaussian distribution describing the transverse profile
131 // is specified by UI command, use value calculated from Moliere Radius
132 if (fSigma < 0) {
133 // 90% of shower is contained within 1 * R_M
134 // 1.645 * std dev of Gaussian contains 90%
135 fSigma = materialRM / 1.645;
136 }
137
138 // Calculate rotation matrix along the particle momentum direction
139 // It will rotate the shower axes to match the incoming particle direction
140 G4RotationMatrix rotMatrix = G4RotationMatrix();
141 double particleTheta = particleDirection.theta();
142 double particlePhi = particleDirection.phi();
143 double epsilon = 1e-3;
144 rotMatrix.rotateY(particleTheta);
145 // do not use (random) phi if x==y==0
146 if (!(std::fabs(particleDirection.x()) < epsilon && std::fabs(particleDirection.y()) < epsilon))
147 rotMatrix.rotateZ(particlePhi);
148
149 // Create hits
150 // First use rejecton sampling to sample from Gamma distribution
151 // then get random numbers from uniform distribution for azimuthal angle, and
152 // from Gaussian for radius
153 G4ThreeVector position;
154 G4double gammaMax = Gamma((fAlpha - 1) / fBeta, fAlpha, fBeta);
155 G4int generatedHits = 0;
156 while (generatedHits < fNbOfHits) {
157 G4double random1 = G4UniformRand() * fLongMaxDepth;
158 G4double random2 = G4UniformRand() * gammaMax;
159 if (Gamma(random1, fAlpha, fBeta) >= random2) {
160 // Generate corresponding rho (phi) from Gaussian (flat) distribution
161 G4double phiPosition = G4UniformRand() * 2 * CLHEP::pi;
162 G4double rhoPosition = G4RandGauss::shoot(0, fSigma);
163 position = particlePosition
164 + rotMatrix
165 * G4ThreeVector(rhoPosition * std::sin(phiPosition),
166 rhoPosition * std::cos(phiPosition), random1 * materialX0);
167 // Create energy deposit in the detector
168 // This will call appropriate sensitive detector class
169 fHitMaker->make(G4FastHit(position, energy / fNbOfHits), aFastTrack);
170 generatedHits++;
171 }
172 }
173 }
174
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176
177 void Par03EMShowerModel::Print() const
178 {
179 G4cout << "Par03EMShowerModel: " << G4endl;
180 G4cout << "Gaussian distribution (transverse plane): \tmu = 0, sigma = "
181 << G4BestUnit(fSigma, "Length") << G4endl;
182 if (fSigma < 0)
183 G4cout << "Negative sigma value means that it will be recalculated "
184 "from the value of the Moliere radius of the detector material, "
185 "taking into account that 90% of the area below the Gaussian "
186 "distribution (from mu - 1.645 sigma to mu + 1.645 sigma) "
187 "corresponds to area within 1 Moliere radius."
188 << G4endl;
189 G4cout << "Gamma distribution (along shower axis): \talpha = " << fAlpha << ", beta = " << fBeta
190 << ", max depth = " << fLongMaxDepth << " X0" << G4endl;
191 if (fAlpha < 0)
192 G4cout << "Negative alpha value means that it will be recalculated "
193 "from the critical energy of the detector material, particle "
194 "type, and beta parameter.\n alpha = beta * T_max, where T_max = "
195 "ln(E/E_C) + C\n where E is particle energy, E_C is critical "
196 "energy in the material, and constant C = -0.5 for electrons and "
197 "0.5 for photons (Eq. (33.36) from PDG)."
198 << G4endl;
199 G4cout << "Number of created energy deposits: " << fNbOfHits << G4endl;
200 }
201