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