| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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(G4Strin
40 : G4VFastSimulationModel(aModelName, aEnvelo
41 fMessenger(new Par03EMShowerMessenger(this
42 fHitMaker(new G4FastSimHitMaker)
43 {}
44
45 //....oooOO0OOooo........oooOO0OOooo........oo
46
47 Par03EMShowerModel::Par03EMShowerModel(G4Strin
48 : G4VFastSimulationModel(aModelName),
49 fMessenger(new Par03EMShowerMessenger(this
50 fHitMaker(new G4FastSimHitMaker)
51 {}
52
53 //....oooOO0OOooo........oooOO0OOooo........oo
54
55 Par03EMShowerModel::~Par03EMShowerModel() = de
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 G4bool Par03EMShowerModel::IsApplicable(const
60 {
61 return &aParticleType == G4Electron::Electro
62 || &aParticleType == G4Positron::Posi
63 || &aParticleType == G4Gamma::GammaDe
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oo
67
68 G4bool Par03EMShowerModel::ModelTrigger(const
69 {
70 // Check energy
71 if (aFastTrack.GetPrimaryTrack()->GetKinetic
72 return false;
73 }
74 // Check length of detector
75 // Calculate depth of the detector along sho
76 // will fit inside. Required max shower dept
77 // can be changed with UI command `/Par03/fa
78 G4double X0 = aFastTrack.GetPrimaryTrack()->
79 auto particleDirection = aFastTrack.GetPrima
80 auto particlePosition = aFastTrack.GetPrimar
81 G4double detectorDepthInMM =
82 aFastTrack.GetEnvelopeSolid()->DistanceToO
83 G4double detectorDepthInX0 = detectorDepthIn
84 // check if detector depth is sufficient to
85 if (detectorDepthInX0 < fLongMaxDepth) {
86 return false;
87 }
88 return true;
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oo
92
93 void Par03EMShowerModel::DoIt(const G4FastTrac
94 {
95 // Remove particle from further processing b
96 aFastStep.KillPrimaryTrack();
97 aFastStep.ProposePrimaryTrackPathLength(0.0)
98 G4double energy = aFastTrack.GetPrimaryTrack
99 // No need to create any deposit, it will be
100 // G4FastSimHitMaker that will call the sens
101 aFastStep.ProposeTotalEnergyDeposited(0);
102 auto particlePosition = aFastTrack.GetPrimar
103 auto particleDirection = aFastTrack.GetPrima
104
105 // Calculate how to create energy deposits
106 // Following PDG 33.5 chapter
107 // material calculation assumes homogeneous
108 auto material = aFastTrack.GetPrimaryTrack()
109 G4double materialX0 = material->GetRadlen();
110 G4double materialZ = material->GetZ();
111 // EC estimation follows PDG fit to solids i
112 G4double materialEc = 610 * MeV / (materialZ
113 // RM estimation follows PDG Eq. (33.37) (rm
114 G4double materialRM = 21.2052 * MeV * materi
115 G4double particleY = energy / materialEc;
116 // Estimate shower maximum and alpha paramet
117 // that describes the longitudinal profile (
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()->GetParti
123 particleTmax += 0.5;
124 }
125 else {
126 particleTmax -= 0.5;
127 }
128 fAlpha = particleTmax * fBeta + 1;
129 }
130 // Unless sigma of Gaussian distribution des
131 // is specified by UI command, use value cal
132 if (fSigma < 0) {
133 // 90% of shower is contained within 1 * R
134 // 1.645 * std dev of Gaussian contains 90
135 fSigma = materialRM / 1.645;
136 }
137
138 // Calculate rotation matrix along the parti
139 // It will rotate the shower axes to match t
140 G4RotationMatrix rotMatrix = G4RotationMatri
141 double particleTheta = particleDirection.the
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()) < eps
147 rotMatrix.rotateZ(particlePhi);
148
149 // Create hits
150 // First use rejecton sampling to sample fro
151 // then get random numbers from uniform dist
152 // from Gaussian for radius
153 G4ThreeVector position;
154 G4double gammaMax = Gamma((fAlpha - 1) / fBe
155 G4int generatedHits = 0;
156 while (generatedHits < fNbOfHits) {
157 G4double random1 = G4UniformRand() * fLong
158 G4double random2 = G4UniformRand() * gamma
159 if (Gamma(random1, fAlpha, fBeta) >= rando
160 // Generate corresponding rho (phi) from
161 G4double phiPosition = G4UniformRand() *
162 G4double rhoPosition = G4RandGauss::shoo
163 position = particlePosition
164 + rotMatrix
165 * G4ThreeVector(rhoPositi
166 rhoPositi
167 // Create energy deposit in the detector
168 // This will call appropriate sensitive
169 fHitMaker->make(G4FastHit(position, ener
170 generatedHits++;
171 }
172 }
173 }
174
175 //....oooOO0OOooo........oooOO0OOooo........oo
176
177 void Par03EMShowerModel::Print() const
178 {
179 G4cout << "Par03EMShowerModel: " << G4endl;
180 G4cout << "Gaussian distribution (transverse
181 << G4BestUnit(fSigma, "Length") << G4
182 if (fSigma < 0)
183 G4cout << "Negative sigma value means that
184 "from the value of the Moliere r
185 "taking into account that 90% of
186 "distribution (from mu - 1.645 s
187 "corresponds to area within 1 Mo
188 << G4endl;
189 G4cout << "Gamma distribution (along shower
190 << ", max depth = " << fLongMaxDepth
191 if (fAlpha < 0)
192 G4cout << "Negative alpha value means that
193 "from the critical energy of the
194 "type, and beta parameter.\n alp
195 "ln(E/E_C) + C\n where E is part
196 "energy in the material, and con
197 "0.5 for photons (Eq. (33.36) fr
198 << G4endl;
199 G4cout << "Number of created energy deposits
200 }
201